# [Math] Runge-Kutta 4 – solving system of 6 differential equations (BVP)

maxima-softwarenumerical methodsordinary differential equationspolynomials

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.$

#### Best Answer

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.