Equazione differenziale difficile da risolvere con Runge-Kutta

[Fede_16]

Salve! Mi sono imbattuto in un'equazione differenziale particolare, probabilmente molto difficile, che vorrei risolvere numericamente con i metodi di Runge-Kutta. Il problema è il seguente, fissato $beta=80$:
\begin{cases}

x'=dx/d\phi=\dfrac{ cos(\phi) }{ 2+\beta \cdot z-(1/x) \cdot sin(\phi) } \\

z'=dz/d\phi=\dfrac{ sin(\phi) }{ 2+\beta \cdot z-(1/x) \cdot sin(\phi) } \\

x(\phi=0)=0; \\ z(\phi=0)=0

\end{cases}

Con $\phi \in [0,\pi/2]$. Non so nemmeno se il problema è ben posto ad essere onesto. Fatto sta che ho provato molto naif con un il metodo di Eulero Esplicito e matlab mi pone errore del codice.

function [t,Y]=euleroesplicitoV2(f,a,b,y0,N)

h=(b-a)/N;
m=length(y0);

t=zeros(N+1,1);
Y=zeros(N+1,m);

t(1)=a;

Y(1,:)=y0;

for i=1:N

t(i+1)=a+i*h;

Y(i+1,:)=Y(i,:) + h*(f(t(i),Y(i,:))); % f deve essere riga
end

Ho dato in input f = @ (phi,x,z) [ (cos(phi))/(2+beta*z-(sin(phi)/x )) ; (sin(phi))/ (2+beta*z-(sin(phi)/x ))], $a=0$, $b=\pi/2$, $y0=[0;0]$ e $N=100$.

Qualcuno ha idea su come muoversi? Grazie per l'attenzione! :]

[sellacollesella]

.

[Fede_16]

Grazie! Su quale programma hai inserito quel codice?
Sull'osservazione per lo zero, adottando la tua soluzione dello "zero ingegneristico", ho provato comunque ad utilizzare il metodo di Runge-Kutta del quarto ordine. Il codice è il seguente, non dà errori che impediscono la risoluzione, ma il risultati che dà forse non sono soddisfacenti.

[bgcolor=lightyellow]function [T,V]=capillaryrisesol(N)

beta=80;

f=@(phi,x,z) [ (cos(phi))/(2+beta*z-(sin(phi)/x) ) ; (sin(phi))/(2+beta*z-(sin(phi)/x ))];

a=0.0000000001;
b=pi/2;

v0=[a ;a];

m=length(v0);

h=( b-0 )/N;

T=zeros(1,N+1);
V=zeros(m,N+1);

V(:,1)=v0;
T(1)=0;

k1=f(T(1),a,a);

k2=f( T(1)+h/2, a + (h/2)*k1(1), a + (h/2)*k1(2) );

k3=f( T(1)+h/2, a + (h/2)*k2(1), a + (h/2)*k2(2) );

k4=f( T(1)+h, a+h*k3(1), a+h*k3(2) );

for n=1:N

T(n+1)=0+n*h;

V(:,n+1)=V(:,n) + (h/6)*(k1+2*k2+2*k3+k4);

end[/bgcolor]

Facendo plot(V(1,:), V(2,:)) esce una retta, ma la soluzione grafica dovrebbe essere effettivamente quella trovata da te. C'è qualcosa che non va nel codice oppure sbaglio nel plottare? Essenzialmente dovrei plottare tutte le coppie date da tutte le colonne di $V$ (che è una matrice 2x1001).

[Fede_16]

UPGRADES: Plottando su scala logaritmica, con

V1=V(1,:); V2=V(2,:);
semilogx(V1,V2)

Esce fuori un grafico che effettivamente ricorda quello trovato. Ma incrementa troppo tardi, con $N=50$, infatti prima di $x=10^2$ circa la soluzione è zero costantemente, mentre dopo inizia a crescere. Dovrebbe crescere molto prima. Probabilmente il $h$ è troppo grande, meglio usare $h^4$. Utilizzando tale passo si ha che la funzione cresce poco prima di $x=10^-3$. Variando $N$ si ottengono risultati diversi.

Sono molto confuso

[sellacollesella]

.

[Fede_16]

Grazie millee! Non sono praticissimo di Mathematica, ma è la volta buona che inizio ad approcciarmi
Comunque, ho ragionato con matlab e sono giunto ad una conclusione che dovrebbe essere equivalente alla tua. Per chiunque sia curioso, lascio qua il codice:

function [T,V]=capillaryrisesol(N)

beta=80;

f=@(phi,x,z) [ (cos(phi))/(2+beta*z-(sin(phi)/x) ) ; (sin(phi))/(2+beta*z-(sin(phi)/x ))];

a=0.000001;
b=pi/2;

v0=[a ;0];

m=length(v0);

T=linspace(0,b,N+1);
T(1)=a;
V=zeros(m,N+1);

V(:,1)=v0;

for n=1:N

    h=T(n+1)-T(n);

    k1=f(T(n),V(1,n),V(2,n));

    k2=f( T(n)+h/2, V(1,n) + (h/2)*k1(1), V(2,n) + (h/2)*k1(2)  );

    k3=f( T(n)+h/2, V(1,n) + (h/2)*k2(1), V(2,n) + (h/2)*k2(2) );

    k4=f( T(n)+h, V(1,n)+h*k3(1),  V(2,n)+h*k3(2) );

    V(:,n+1)=V(:,n) + (h/6)*(k1+2*k2+2*k3+k4);

end

Con

[T,V]=capillaryrisesolV2(50);
V1=V(1,:); V2=V(2,:);
plot(V1,V2)

Ottengo:



Che mi pare sia piuttosto simile alla tua soluzione

[sellacollesella]

.

[Fede_16]

Perfetto sono equivalenti! Chiaramente qualche cifra decimale cambia però con scarto piuttosto basso. Immagino sia dovuto alla differenza tra i due metodi :] Ho verificato anche sul libro e i risultati son proprio questi Grazie ancora!

[feddy]

Ho visto questo post solo ora. Formalmente, una singolarità presente nel problema continuo non si può eliminare impunemente nel discreto. L'approccio di usare una quantità vicino all'epsilon macchina come $0$ funziona perché la singolarità è del tipo $\frac{\sin(x)}{x}$, e per $x \rightarrow 0$ il limite è finito (e notevole ).

A proposito, se provate a prendere uno schema implicito, ad esempio Eulero implicito si ha la seguente ricorrenza:

$$[x_{k+1},z_{k+1}] = [x_{k},z_{k}] + hF([x_{k+1},z_{k+1}])$$ dove $h$ è il passo temporale, preso uniforme per comodità, e $F$ il membro di destra della ODE. La comodità di questo schema è che la valutazione di $F$ non viene fatta in $\phi=0$ direttamente (ma va risolto ad ogni istante temporale un sistema lineare).

[Fede_16]

Grazie per l'approfondimento! :]] Se non ci fosse stato il fortuito limite finito, il problema risultava mal posto e quindi non risolubile?

[feddy]

Uno schema numerico per ODE ha bisogno di alcune ipotesi che valgano nel problema continuo. Ad esempio, la convergenza, etc si dimostra espandendo in serie di Taylor, etc... Se questo cade, non c'e' schema numerico che possa sanare la situazione in questo contesto.

[Fede_16]

Capisco, grazie molte! :]