Problemas con las derivadas usando Newton-Raphson en MatLab

Question

Me está resultando muy difícil entender cuál es la mejor manera de expresar el siguiente sistema de ecuaciones en MatLab para poder resolverlo. Las ecuaciones provienen de la solución de similitud de Von Karman para Navier-Stokes para un flujo de disco giratorio.
He tomado las ecuaciones originales y las he convertido en un conjunto de cinco EDOs de primer orden, pero no consigo averiguar cómo traduzco las ecuaciones que tengo a un formato que funcione en MatLab, dado que todos los ejemplos que he visto vienen en la forma $a x^2 + b x + c = 0$ , o algún equivalente, no sé los diferenciales entran en esto: ¿sólo tomo el lado derecho de la ecuación, o los traigo y lo hago todo igual a cero, o qué?

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

Con condiciones de contorno $U(0) = 0, V(0) = 0, W(0) = 0, U(20) = 0, V(20) = -1 $
y las estimaciones iniciales para $g_1(0) = 0.52$ y $g_2(0) = -0.61$

Las funciones U, V y W se introducen en las ecuaciones siguientes para producir los tres componentes de la velocidad del flujo:

(6) $u = r\Omega *U(n)$ ,
(7) $v= r\Omega *V(n)$ ,
(8) $v = \sqrt{ \nu\Omega} * W(n)$
donde $n = z\sqrt\frac{\nu} {\Omega}$ y $\nu=$ viscosidad, $\Omega =$ velocidad angular, $z =$ la coordenada axial y $r =$ coordenada radial

El segundo conjunto de ecuaciones, (6),(7) y (8) lo he añadido simplemente por claridad, para que se pueda entender qué es lo que realmente estoy buscando.
Así que sí- pregunta principal: ¿Cómo escribo realmente las ecuaciones 1-5 en un formato que se ajuste a un método de rodaje en MatLab?

Muchas gracias, Mike

Answer 1

Por favor, no abra varios mensajes para la misma pregunta. Trate de entender el respuestas dado y pedir, si es necesario, más información en los comentarios.

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

A continuación, haz una función que conecte los valores iniciales con dos parámetros con los valores finales

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

y luego aplicas algún tipo de método Newton (probado con scilab, comprueba la sintaxis correcta de matlab)

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

Entonces la iteración termina estacionaria en u1=0.5106376 y v1=-0.6089510 . Sin embargo, los valores iniciales que conducen a un diagrama como en el papel son u1=0.5102326 y v1=-0.6159221 , siguen siendo válidos también para intervalos de integración más largos, al contrario que el primer par. Esta es la sensibilidad mencionada en el documento.

Hay que tener cuidado con las velocidades iniciales, el sistema sólo tiene un rango muy estrecho en el que no explota antes de alcanzar $t=5$ .