Solving a system of two second order ODEs using Runge-Kutta method (ode45) in MATLAB

Question

I'm trying to solve the system of differential equations outlined in Von Karman's rotating disk flow. I got them into a system of ordinary differential equations:

F(n), G(n), H(n)

$$F'' = -G^2 + F^2 - F'H$$ $$G'' = 2FG - HG'$$ $$H' = -2F$$

with boundary conditions $F(0) = 0$, $G(0) = 1$, $H(0)=1$, $F(\infty) = 0$, $G(\infty) = 0$.

And then converted that into a system of 1st order ODEs featuring $F$, $F_2$, $G$ and $G_2$:

$$\begin{align*}F' &= F_2\\ G' &= G_2\\ F_2' &= -G^2 + F^2 - F_2F^2\\ G_2' &= 2FG - F^2G_2 \end{align*}$$

However here is where I've hit a wall. I first tried solving them analytically using DSolve, however that failed to return an explicit solution, so after some reading I discovered that I could try using a numerical shooting method. Google seems to indicate that ode45 is my best bet, however I'm not particularly adept at MatLab, and I've been tearing my hair out trying to understand how to implement it successfully.

From what I've read in other papers, some have had success using initial conditions of F2'(0) = 0.52 and G2'(0) = -0.61 and a replacement of infinity with 20, if that's any help at all.

Bottom line is, I really just can't seem to get my head around how to solve this in MatLab. If anyone could show me how to, or at least demonstrate the proper syntax for describing this type of system, I'd hugely appreciate it.

Thanks for reading

Mike

Answer 1

Are you sure about your elimination of $H$? Could you explain how you come to replace $H=-F^2$, because that implies $H'=-2F·F'$, which is still quadratic in $F$. Quasi-archeological reconstruction suggests that the following procedure is closest to the original system in the paper.

function doty = odefunc(t,y)
    F=y(1); F1=y(2); G=y(3); G1=y(4); H=y(5)
    doty = [ F1, -G^2 + F^2 - F1*H,
             G1, 2*F*G - H*G1,
             2*F ]
end

This now is a dimension 5 system, thus requires 5 initial values.

---

Using the structure of the system in the comments, and the other question closer to the original paper, one can reconstruct that the correct system has indeed the structure $$\begin{align} Z'' &= Z^2-HZ',&~~~ Z=F+iG&=U+i(V+1), \\ H'&=2F,& ~~~ H&=-W \end{align}$$ as one could expect from a rotational setup.

Using a boundary value solver this gives a solution over the interval $[0,8]$ that has first and final vector

[-7.47756897e-23  1.00000000e+00  5.10116328e-01 -6.15849164e-01  0.00000000e+00]
[ 0.00000000e+00  0.00000000e+00 -7.04991619e-04 -9.34891188e-04  8.74618873e-01]

with a plot

However, the python script fails fatally for larger intervals, so convergence towards a solution with right boundary at infinity is unclear.