Trouble with derivatives using Newton-Raphson in MatLab

Question

I'm finding it very difficult to get my head around how best to express the following system of equations in MatLab in order to solve it. The equations come from Von Karman's similarity solution to Navier-Stokes for a rotating disk flow.
I've taken the original equations and turned them into a set of five 1st order ODEs, but I can't work out how I translate the equations I have into a format that works in MatLab, given that all the examples I've seen come in the form $a x^2 + b x + c = 0$, or some equivalent, I don't know the differentials come into this: am I only taking the right hand side of the equation, or do I bring them over and make it all equal to zero, or what?

(1) $g_1 = U'$

(2) $g_2 = V'$

(3) $W' = -2U$

(4) $g_1' = U^2 - (V+1)^2 + g_1 W$

(5) $g_2' = 2U(V+1)+g_2 W$

With boundary conditions $U(0) = 0, V(0) = 0, W(0) = 0, U(20) = 0, V(20) = -1 $
and initial guesses for $g_1(0) = 0.52$ and $g_2(0) = -0.61$

The functions U, V, and W themselves then feed into the equations below in order to produce the three components of velocity for the flow:

(6) $u = r\Omega *U(n)$,
(7) $v= r\Omega *V(n)$,
(8) $v = \sqrt{ \nu\Omega} * W(n)$
where $n = z\sqrt\frac{\nu} {\Omega}$ and $\nu=$ viscosity, $\Omega =$ angular velocity, $z =$ the axial coordinate and $r =$ radial coordinate

The second set of equations, (6),(7) and (8) I have simply added for clarity, so you can understand what it is I'm actually looking for.
So yes- main question: How do I actually write equations 1-5 in a format that works fits a shooting method in MatLab?

Many thanks, Mike

Answer 1

Please do not open multiple posts for the same question. Try to understand the answers given and ask, if needed, for more information in the comments.

function doty = odefunc(t,y)
    U=y(1); U1=y(2); V=y(3); V1=y(4); W=y(5)
    doty = [ U1; U^2 - (V+1)^2 + U1*W;
             V1; 2*U*(V+1) + W*V1;
             -2*U ]
end

Then make a function that connects the initial values with two parameters to the end values

function uv20 = F(u10,v10)
    t0 = 0; tf = 20;
    [t, y] = ode45(odefunc, [t0 tf], [ 0; u10; 0; v10; 0])
    uv20 = [ y(1,2); y(3,2) ]
end

and then you apply some kind of Newton method (tested with scilab, check for correct matlab syntax)

u1=0.52; v1=-0.61

h=1e-5;

for k=1:10 do
    Fval = F(u1,v1)
    dFdu1 = (F(u1+h, v1) - F(u1-h, v1))/(2*h);
    dFdv1 = (F(u1, v1+h) - F(u1, v1-h))/(2*h);

    dF = [ dFdu1 dFdv1 ] 

    uvnext = [ u1;v1 ] - dF^(-1)*(Fval-[0;-1])

    u1 = uvnext(1); v1 = uvnext(2);

end

Then the iteration ends up stationary at u1=0.5106376 and v1=-0.6089510. However, the initial values leading to a diagram like in the paper are u1=0.5102326 and v1=-0.6159221, they remain valid also for longer integration intervals, contrary to the first pair. This is the sensitivity mentioned in the paper.

Be careful with the initial velocities, the system only has a very narrow range where it does not explode before reaching $t=5$.