6

I'm facing a tricky problem. I need to solve a system of 6 differential equations numerically, but I don't have 6 IVP (initial value problem) conditions, instead I have 6 BVP (boundary valye problem) conditions (3 conditions at $x=0$ and 3 conditions at $x=l$, being $l$ the last point to compute, thus $0<x<l$).

So I need to convert the 3 boundary conditions at $x=l$ into 3 at $x=0$ (through the shooting method for example) so I have all 6 IVPs to solve my problem numerically through RK4. So far so good.

I have two problems to solve. The first has the following form:

\begin{cases} \frac{dT}{dx}&=f_1(g) &(1)\\ \frac{dV}{dx}&=f_2(e) &(2)\\ \frac{dM}{dx}&=f_3(V,g) &(3)\\ \frac{dK}{dx}&=f_4(x,M,T) &(4)\\ \frac{dg}{dx}&=f_5(x,T) &(5)\\ \frac{de}{dx}&=f_6(K) &(6)\\ \end{cases}

With the following BVPs:

\begin{cases} T(0)&=c_1 \\ T(l)&=0 \\ V(0)&=c_2 \\ V(l)&=0 \\ M(0)&=c_3 \\ M(l)&=0 \\ \end{cases}

Through the shooting method I need to find $K(0)$, $g(0)$ and $e(0)$. In this case I can trivially find $g(0)$ because (1) and (5) are linearly dependent. To find $K(0)$ and $e(0)$ I shoot, for example:

\begin{cases} K(0)=+1 \text{ and } e(0)=+1 \\ K(0)=+1 \text{ and } e(0)=-1 \\ K(0)=-1 \text{ and } e(0)=-1 \\ K(0)=-1 \text{ and } e(0)=+1 \\ \end{cases}

Then I record the final values for $V(l)$ and $M(l)$ for each case, make an approximating polynomial for each ($f_V=f(K,e)$ and $f_M=f(K,e)$ respectively), and solve the resulting system for $K$ and $e$:

\begin{cases} f_V(K,e)=0 \\ f_M(K,e)=0 \end{cases}

The resulting $K(0)$ and $e(0)$ give me fairly good values (if I use approximating polynomials of order 1 I have an error of $10^-8$, more than acceptable for what I need). Thus, my first problem is solved.

The second problem is the real issue. The system is now:

\begin{cases} \frac{dT}{dx}&=f_1(e,g)\\ \frac{dV}{dx}&=f_2(e,g) \\ \frac{dM}{dx}&=f_3(V,e,g) \\ \frac{dK}{dx}&=f_4(x,M,T) \\ \frac{dg}{dx}&=f_5(x,T) \\ \frac{de}{dx}&=f_6(K) \\ \end{cases}

Also, there are now hyperbolic tangents involved (I have $g(x)$ inside $tanh()$ for example). I can't find any of $K(0)$, $g(0)$ or $e(0)$ trivially because now all equations depend on the others. So, my initial idea was to use the same shooting methodology but creating polynomials with 3 variables ($K$, $g$ and $e$) instead of just $K$ and $e$ like I used before. This isn't producing viable results, even if I use higher orders of polynomials.

My question is: what can I do ? I don't mind changing methodologies or solving using other methods, I just need some directions so I can correctly solve this. I wholeheartedly welcome different perspectives on this. I'm using maxima, for what it's worth.

Also, I can determine the correct $K(0)$, $g(0)$ and $e(0)$ (through an external Maple file), and if I input those into my program I get correct results, so I know all the equations are correctly inputted.

Thank you !

EDIT2:

The system I'm aiming to solve is the following:

\begin{cases} \frac{dT(x)}{dx}&={{3339\,{\it g(x)}\,\tanh \left({{50000\,\sqrt{8427\,{\it g(x)}^2+30036 \,{\it e(x)}^2}}\over{163611}}\right)}\over{\sqrt{8427\,{\it g(x)}^2+ 30036\,{\it e(x)}^2}}}\\ \frac{dV(x)}{dx}&={{12600\,{\it e(x)}\,\tanh \left({{50000\,\sqrt{8427\,{\it g(x)}^2+ 30036\,{\it e(x)}^2}}\over{163611}}\right)}\over{\sqrt{8427\,{\it g(x)}^ 2+30036\,{\it e(x)}^2}}} \\ \frac{dM(x)}{dx}&={{5\,\sqrt{8427\,{\it g(x)}^2+30036\,{\it e(x)}^2}\,{\it V(x)}-10017\, {\it g(x)}\,\tanh \left({{50000\,\sqrt{8427\,{\it g(x)}^2+30036\, {\it e(x)}^2}}\over{163611}}\right)}\over{5\,\sqrt{8427\,{\it g(x)}^2+ 30036\,{\it e(x)}^2}}} \\ \frac{dK(x)}{dx}&=-{{60092431565\,x+15000000000\,{\it T(x)}+25000000000\,{\it M(x)}- 2410745228515}\over{48076923076923}} \\ \frac{dg(x)}{dx}&={{3281046763449\,x+1069000000000\,{\it T(x)}-156626689476919}\over{ 5250000000000000}} \\ \frac{de(x)}{dx}&={\it K(x)} \\ \end{cases}

Boundary conditions are:

\begin{cases} T(0)=200 \\ V(0)=-4.8073945252 \\ M(0)=-47.1403817188 \\ T(l)=0 \\ V(l)=0 \\ M(l)=0 \\ \end{cases}

Where $l=25mm.$

Marcelo
  • 153
  • Since you mentioned Maple, you might note that Maple's dsolve includes solvers for boundary value problems. Unless there's some particular reason you have to use RK4, you're probably better off giving the problem to Maple as a boundary value problem. – Robert Israel Mar 15 '13 at 23:50
  • If you have BVP why don't you use finite difference method? – tom Mar 15 '13 at 23:56
  • Do you have different functions? Otherwise I'm not sure I know how one function can accept different number of arguments. – Kaster Mar 15 '13 at 23:58
  • @RobertIsrael I can use Maple for verification purposes, but use maxima for the actual calculations as I mentioned. That is necessary. – Marcelo Mar 16 '13 at 08:26
  • @tom I still have the same problem of not having $K(0)$, $g(0)$ and $e(0)$. Without them I can't start the finite difference analysis on those equations, I think ? – Marcelo Mar 16 '13 at 08:32
  • @Kaster I don't think I understand what you mean, but each of the 6 functions in the system are different yes. – Marcelo Mar 16 '13 at 08:34
  • @Marcelo it's just confusing that you denote all of them as $f$. And it's hard to see that (1) and (5) are linearly dependent. – Kaster Mar 16 '13 at 19:48
  • @Kaster you are absolutely right, I didn't notice that until you pointed it out. I edited the post to clarify that point. – Marcelo Mar 16 '13 at 20:38

1 Answers1

1

Ok I try to answer but I have no clue if it works. I'll make simpler case where you solve:

$$ u' = f(u,v,x) \\ v'=g(u,v,x)$$

But you only know u(0) and u(l).

Denote grid points $0=x_0 < x_1 < \dots < x_n = l$ with $x_{i+1}-x_i = h$. $u(x_i) = u_i$.

Now using central difference we arrive to system of equations. $$ \frac{1}{2h} \begin{pmatrix} 0 & 1 & 0 & \cdots & 0&0 \\ -1& 0 & 1 & \cdots & 0&0 \\ 0 & -1 & 0 & \cdots & 0&0 \\ \vdots & \vdots &\vdots& \ddots & \vdots &\vdots \\ 0& 0&0& \cdots &0& 1\\ 0&0&0&\cdots &-1 & 0 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\u_3 \\ \vdots \\u_{n-2} \\ u_{n-1} \end{pmatrix} + \frac{1}{2h} \begin{pmatrix} -u(0) \\ 0 \\ 0 \\ \vdots \\0 \\ u(l) \end{pmatrix} =\begin{pmatrix} f(u_1,v_1,x_1) \\ f(u_2,v_2,x_2) \\ f(u_3,v_3,x_3) \\ \vdots \\f(u_{n-2},v_{n-2},x_{n-2}) \\ f(u_{n-1},v_{n-1},x_{n-1}) \end{pmatrix} $$

$$ \frac{1}{2h} \begin{pmatrix} -2 & 2 & 0 & \cdots & 0&0 \\ -1& 0 & 1 & \cdots & 0&0 \\ 0 & -1 & 0 & \cdots & 0&0 \\ \vdots & \vdots &\vdots& \ddots & \vdots &\vdots \\ 0& 0&0& \cdots &0& 1\\ 0&0&0&\cdots &-2 & 2 \end{pmatrix} \begin{pmatrix} v_0 \\ v_1 \\v_3 \\ \vdots \\v_{n-1} \\ v_{n} \end{pmatrix} =\begin{pmatrix} g(u_0,v_0,x_0) \\ g(u_1,v_1,x_1) \\ g(u_2,v_2,x_2) \\ \vdots \\g(u_{n-1},v_{n-1},x_{n-1}) \\ g(u_{n},v_{n},x_{n}) \end{pmatrix} $$

When discretizing second equation I used forward and backward difference at point 0 and l respectively.

Now you have system of nonlinear equations which you can try to solve.

tom
  • 4,737
  • 1
  • 23
  • 43
  • This helps if I try to solve my first problem, but unfortunately it doesn't help with my second problem. The resulting system is too complex (even if I use just 4 steps, so $h=l/4$) to solve. The reason must be the existence of functions $e(x)$ and $g(x)$ inside a square root and a hyperbolic tangent in the first 3 equations of my system. – Marcelo Mar 16 '13 at 18:10
  • What do you mean byt "too complex"? and sorry but I don't understand your reasoning "The reason must be ..." at all. – tom Mar 16 '13 at 20:51
  • And write down exact form of $f_1,\cdot,f_6$. I would like to try to solve that system of equations. – tom Mar 16 '13 at 20:54
  • by too complex I meant both maxima and Maple errored out without solutions after some time. I updated the main post with the complete system I'm trying to solve ! – Marcelo Mar 16 '13 at 23:03
  • Any luck Tom ? Thanks. – Marcelo Mar 18 '13 at 13:54
  • sorry I didn't have time yet – tom Mar 19 '13 at 23:21
  • It's OK, I appreciate your help when you can ! Thanks again. – Marcelo Mar 19 '13 at 23:32