Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi

Messaggioda Zurtolo » 17/07/2019, 15:53

Buongiorno a tutti, dovrei risolvere numericamente un sistema di equazioni differenziali dove la prima equazione è del terzo ordine e la seconda equazione è del secondo ordine.
Il sistema è il seguente

$ w'''+1/rw''-1/r^2w'=12/h^2w'*( u'+vu/r+1/2*(w')^2 )+(qr)/(2*D) $
$ u''+1/ru'-u/r=-((1-v)/(2r))(w')^2-w'w'' $

Dove:
* r è la variabile dipendente
* h=1 [millimetri]
* E=200000 [MegaPascal]
* v=0.3
* $ D=(Eh^3)/(12(1-v^2)) $
* q=20 [MegaPascal]

(Il sistema fisicamente rappresenta la deflessione di una piastra circolare di acciaio di raggio a e spessore h soggetta a un carico uniformemente distribuito q; w e u sarebbero rispettivamente la deflessione e lo spostamento radiale della piastra in funzione della posizione radiale r )

L'integrazione deve avvenire sull'intervallo [a b], dove:
*a=12
*b=0

Con le seguenti condizioni al contorno:
*w(r=a)=0
*w'(r=a)=0
*w'(r=0)=0
*u(r=a)=0
*u'(r=0)=0
(Tali condizioni fisicamente rappresentano un'incastro.)

Purtroppo le condizioni al contorno non sono definite in un unico punto, ma in due. Per risolvere la cosa ho pensato di utilizzare il metodo di shooting trasformando il problema in un sistema di ODE. Ho quindi prima testato con successo un codice di calcolo in linguaggio Matlab per una singola equazione (in particolare ho risolto un equazione di una discussione già riportata nel forum (https://www.matematicamente.it/forum/vi ... p?t=189081) e poi adattato il codice per questo sistema. Il codice funziona e calcola correttamente (ho confrontato con analisi FEM) ma non se h diventa troppo piccolo (per esempio h=0.1) o se q diventa troppo grande perchè la soluzione oscilla. Mentre invece usando la funzione di matlab per i problemi con i valori ai bordi "bvp4c" si ottiene la soluzione.
Volevo implementare un'altro schema (ad esempio uno schema di shooting multiplo) per ovviare a questo inconveniente. Ho provato a implementare uno shooting multiplo ma non ho ben capito come costruirlo. Chiedo quindi gentilmente se qualcuno può aiutarmi. Allego anche gli script che ho steso. Grazie a tutti un saluto

Codice:
function shooting_method
clc
clear all
close all

global a h E v q D

a=12;
h=0.1;
E=200000;
v=0.3;
q=20;

D=E*h^3/12/(1-v^2);

x=[0 0];
options=optimset('Display','iter');
x1=fsolve(@solver,x)

end

function F=solver(x)
global a h E v q D

%options=odeset('RelTol', 1e-08, 'AbsTol', [1e-08 1e-08 1e-08 1e-08 1e-08]);
[t,u]=ode(@equation, [a 0.01], [0 0 x(1) 0 x(2)]);

s=length(t);
F=[ abs( u(s,2) ) , abs( u(s,4) )  ]

clf
plot(t,u(:,1),'-or')
xlabel(' r[mm] ');
ylabel(' dy [mm] ');

end


function dy=equation(r,y)
global a h E v q D
dy=zeros(5,1);

dy(1)=y(2);
dy(2)=y(3);
dy(3)=-y(3)/r+y(2)/r^2+12/h^2*y(2)*( y(5)+v*y(4)/r+0.5*y(2)^2 )+q*r/2/D;
dy(4)=y(5);
dy(5)=-y(5)/r+y(4)/r^2-(1-v)/2/r*y(2)^2-y(2)*y(3);

end


Codice per (https://www.matematicamente.it/forum/vi ... p?t=189081)
Codice:
function shooting_method
clc
clear all
x=0.5;
display('La condizione iniziale equivalente vale')
x1=fzero(@solver,x)

end

function F=solver(x)
options=odeset('RelTol', 1e-08, 'AbsTol', [1e-08, 1e-08]);
a=1;
b=3;
y0=[17 x];
[t,u]=ode45(@equation, [a b], y0, options);
s=length(t);
F=u(s,1)-43/3;
figure(1)
clf
plot(t,u(:,1),'-*r') % plot della soluzione numerica
hold on
plot(t,(t.^2+16./t),'-b') % plot della soluzione esatta

end


function dy=equation(t,y)

dy=zeros(2,1);
dy(1)=y(2);
dy(2)=1/8.*( 32+2.*t.^3-y(1).*y(2) );


end
Zurtolo
Starting Member
Starting Member
 
Messaggio: 1 di 14
Iscritto il: 21/06/2019, 15:32

Re: Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi

Messaggioda Raptorista » 18/07/2019, 10:17

Ciao, un paio di precisazioni innanzitutto

Zurtolo ha scritto:* r è la variabile dipendente

Indipendente?
Zurtolo ha scritto:L'integrazione deve avvenire sull'intervallo [a b], dove:
*a=12
*b=0

\(b < a\)? :?
Zurtolo ha scritto:Con le seguenti condizioni al contorno:
*w(r=a)=0
*w'(r=a)=0
*w'(r=0)=0
*u(r=a)=0
*u'(r=0)=0
(Tali condizioni fisicamente rappresentano un'incastro.)

In questo caso, \(r = 0\) sarebbe \(r = b\)?

C'è un motivo per cui vuoi usare il metodo di shooting? Forse col metodo delle differenze finite si fa prima.

I problemi di oscillazione sono probabilmente solo spazzatura numerica, che matlab ripulisce con qualche trucco dietro le quinte.
Un matematico ha scritto:... come mia nonna che vuole da anni il sistema per vincere al lotto e crede che io, in quanto matematico, sia fallito perché non glielo trovo


Immagine
Avatar utente
Raptorista
Moderatore
Moderatore
 
Messaggio: 5287 di 9616
Iscritto il: 28/09/2008, 19:58

Re: Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi

Messaggioda Zurtolo » 18/07/2019, 16:52

Grazie per l'aiuto.

r sarebbe la variabile indipendente scusa hai ragione. Ho sbagliato a scrivere.

Inoltre non ho precisato nel precedente post, che nel codice ho ricondotto le due equazioni (di cui la prima di grado 3 e la seconda di grado 2) in un sistema di 5 equazioni di grado 1 prima sostituendo:
$ y(1)=w $
$ y(2)=w' $
$ y(3)=w'' $
$ y(4)=N $
$ y(5)=N' $
e poi esplicitando
$ y'(1)=y(2) $
$ y'(2)=y(3) $
$ y'(3)=-(y(3))/r+(y(2))/r^2+12/h^2*y(2)*( y(5)+v*(y(4))/r+0.5*y(2)^2 )+(q*r)/(2*D) $
$ y'(4)=y(5) $
$ y'(5)=-(y(5))/r+(y(4))/r^2-(1-v)/(2r)*y(2)^2-y(2)*y(3) $

Allora b<a nel senso che nel mio codice inizio integro la soluzione al contrario (parto da "a" che sarebbe il raggio esterno della piastra e arrivo a b che è il centro della piastra dove ho la condizione di assialsimmetria). Volendo si può fare anche al contrario. Si r=0 sarebbe quindi r=b.

Non c'è un motivo per cui voglio usare il metodo di shooting se non che essendo un sistema di equazioni differenziali potevo appoggiarmi ai solutori ODE. Qualsiasi metodo che mi permetta di scrivere un solutore "mio" è ben accetto. Sul libro di Calcolo c'è scritto come fare nel caso di una equazione ma nel caso di un sistema come mi devo muovere? Dovrei prima ricondurmi ad un sistema con due equazioni dello stesso grado?

Comunque per ovviare alle oscillazioni ho provato anche solutori per problemi "stiff" come ODE15s ma non riesco a ottenere buoni risultati.

In attesa di risposta vi ringrazio
Zurtolo
Starting Member
Starting Member
 
Messaggio: 2 di 14
Iscritto il: 21/06/2019, 15:32

Re: Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi

Messaggioda Zurtolo » 19/07/2019, 17:31

Per prendere la strada delle differenze finite come devo impostare le equazioni? Ho chiaro come si fa nel caso di una singola equazione ma per un sistema? Devo scalarle allo stesso ordine e scriverle in forma vettoriale?
Zurtolo
Starting Member
Starting Member
 
Messaggio: 3 di 14
Iscritto il: 21/06/2019, 15:32

Re: Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi

Messaggioda Raptorista » 19/07/2019, 17:57

Zurtolo ha scritto: $ y'(1)=y(2) $
$ y'(2)=y(3) $
$ y'(3)=-(y(3))/r+(y(2))/r^2+12/h^2*y(2)*( y(5)+v*(y(4))/r+0.5*y(2)^2 )+(q*r)/(2*D) $
$ y'(4)=y(5) $
$ y'(5)=-(y(5))/r+(y(4))/r^2-(1-v)/(2r)*y(2)^2-y(2)*y(3) $

Non ho controllato i conti, ma concettualmente dovrebbe essere corretto.
Zurtolo ha scritto:Non c'è un motivo per cui voglio usare il metodo di shooting se non che essendo un sistema di equazioni differenziali potevo appoggiarmi ai solutori ODE. Qualsiasi metodo che mi permetta di scrivere un solutore "mio" è ben accetto. Sul libro di Calcolo c'è scritto come fare nel caso di una equazione ma nel caso di un sistema come mi devo muovere? Dovrei prima ricondurmi ad un sistema con due equazioni dello stesso grado?

Premesso che non ho mai usato questo metodo, a quanto ho capito l'idea è di usare un metodo di "root-finding" per trovare la giusta condizione iniziale sulle derivate che ti fa uscire la condizione che ti serve al bordo opposto.
In questo caso hai scambiato la condizione \(w'(b) = 0\) con \(w''(a) = c_1\) e \(u'(b) = 0\) con \(u'(a) = c_2\), quindi devi andare a caccia di \(c_1\) e \(c_2\) tali che la soluzione corrispondente abbia \(w'(b) = 0\) e \(u'(b) = 0\).

Zurtolo ha scritto:Per prendere la strada delle differenze finite come devo impostare le equazioni? Ho chiaro come si fa nel caso di una singola equazione ma per un sistema? Devo scalarle allo stesso ordine e scriverle in forma vettoriale?

Per un sistema hai solo più gradi di libertà. Non serve "scalarle", qualunque cosa significhi. Semplicemente discretizza le derivate con dei rapporti incrementali appropriati.
Un matematico ha scritto:... come mia nonna che vuole da anni il sistema per vincere al lotto e crede che io, in quanto matematico, sia fallito perché non glielo trovo


Immagine
Avatar utente
Raptorista
Moderatore
Moderatore
 
Messaggio: 5292 di 9616
Iscritto il: 28/09/2008, 19:58

Re: Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi

Messaggioda Zurtolo » 22/07/2019, 08:31

Buongiorno, ho sostituito ma non ho ben chiaro come faccio a costruire la matrice per risolvere con le differenze finite il sistema.
Per prima cosa ho posto: $ z=w' $ così avrò un sistema di due equazioni del secondo ordine. w lo calcolerò con un integrazione successiva se riuscirò a venirne fuori.
Utilizzando la differenza finita centrale ottengo:

1) $ (z^(i+1)-2z^(i)+z^(i-1))/(p^2)-(z^(i+1)-z^(i-1))/(2p)*1/r^i-z^i/(r^i)^2=12/h^2*(z^(i+1)-z^(i-1))/(2p)*[ (u^(i+1)-u^(i-1))/(2p)+v*u^i/r^i+1/2*(z^i)^2 ]+qr^i/(2D) $

2) $ (u^(i+1)-2u^(i)+u^(i-1))/(p^2)-(u^(i+1)-u^(i-1))/(2p)*1/r^i-u^i/(r^i)^2=-(1-v)/(2*r^i)*(z^i)^2-z^i*(z^(i+1)-z^(i-1))/(2p) $

dove p è il passo della differenza finita e l'indice i rappresenta l'i-esimo nodo. Non ho capito come faccio a accoppiare queste due formule con un unica matrice. Ho notato che il primo membro di entrambe le equazioni si può associare alla classica matrice tridiagonale, ma il resto come lo sistemo??

Grazie mille per l'aiuto un saluto
Zurtolo
Starting Member
Starting Member
 
Messaggio: 4 di 14
Iscritto il: 21/06/2019, 15:32

Re: Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi

Messaggioda Raptorista » 22/07/2019, 10:45

Bene, sei sulla strada giusta!
Adesso, ipotizzando che tu abbia \(N\) punti nell'intervallo \([a,b]\), avrai un valore \(z^i\) ed un valore \(u^i\) per ciascuno di questi punti, quindi devi costruire un sistema di ordine \(2N\) in cui, per esempio, le prime \(N\) righe contengono le equazioni per le variabili \(z^i\) e le ultime \(N\) righe contengono le equazioni per le variabili \(u^i\). Nel frattempo devi anche assemblare il termine noto, anch'esso di dimensione \(2N\).
Per i nodi di bordo non puoi scrivere un'equazione partendo da quello schema perché ti mancano i nodi a fianco, quindi lì devi imporre le condizioni al bordo come equazione.

In questo caso, e non me ne ero accorto fino ad ora [:-D] il sistema è nonlineare, il che significa che non avrai una matrice ma una generica funzione nonlineare, che poi puoi dare in pasto a una funzione matlab che implementa Newton o qualcosa del genere.
Il fatto che il sistema sia nonlineare rende la cosa leggermente più articolata, ma dovresti farcela senza troppi problemi. Secondo me rimane più semplice del metodo di shooting.
Un matematico ha scritto:... come mia nonna che vuole da anni il sistema per vincere al lotto e crede che io, in quanto matematico, sia fallito perché non glielo trovo


Immagine
Avatar utente
Raptorista
Moderatore
Moderatore
 
Messaggio: 5294 di 9616
Iscritto il: 28/09/2008, 19:58

Re: Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi

Messaggioda Zurtolo » 22/07/2019, 13:57

Ciao grazie mille per la risposta.
Ti chiedo gentilmente se puoi controllare se quello che ho capito è corretto. Io ho capito che:

*Non devo creare alcuna matrice come diciamo nel caso classico delle differenze finite (es. problema di Poisson)

**Per quantor riguarda i nodi "interni". Devo trovare la soluzione per il nodo i-esimo. Ho disponibili le soluzioni al nodo i-1 e i-2 e risolvendo il sistema (non lineare) formato dalle formule "incrementali" che ho trovato nella mia precedente risposta, trovo i valori di u e di z al nodo i.

***Per i valori al bordo non ho capito invece come devo fare...

Ti ringrazio per il tuo tempo, un saluto
Zurtolo
Starting Member
Starting Member
 
Messaggio: 5 di 14
Iscritto il: 21/06/2019, 15:32

Re: Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi

Messaggioda Raptorista » 22/07/2019, 14:03

Zurtolo ha scritto:*Non devo creare alcuna matrice come diciamo nel caso classico delle differenze finite (es. problema di Poisson)

Non devi creare una matrice perché l'equazione è nonlineare. Se l'equazione fosse lineare potresti creare un sistema lineare come si fa di solito.
Zurtolo ha scritto:**Per quantor riguarda i nodi "interni". Devo trovare la soluzione per il nodo i-esimo. Ho disponibili le soluzioni al nodo i-1 e i-2 e risolvendo il sistema (non lineare) formato dalle formule "incrementali" che ho trovato nella mia precedente risposta, trovo i valori di u e di z al nodo i.

Non hai \(u^{i-1}\), tutti i valori sono incogniti e li devi trovare tutti insieme risolvendo il sistema nonlineare.
Zurtolo ha scritto:***Per i valori al bordo non ho capito invece come devo fare...

Se la condizione al bordo dice \(u(b) = 5\) allora per il valore \(u^N\), che suppongo essere quello situato in \(x = b\), dovrai scrivere l'equazione \(u^N = 5\). Similmente per gli altri.
Un matematico ha scritto:... come mia nonna che vuole da anni il sistema per vincere al lotto e crede che io, in quanto matematico, sia fallito perché non glielo trovo


Immagine
Avatar utente
Raptorista
Moderatore
Moderatore
 
Messaggio: 5297 di 9616
Iscritto il: 28/09/2008, 19:58

Re: Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi

Messaggioda Zurtolo » 25/07/2019, 12:38

Buongiorno, ti ringrazio per la risposta. Allora forse ho capito , ti chiedo ancora una volta gentilmente di verificare se il procedimento che ho capito va bene...
Ipotizziamo per semplificità di discretizzare con 4 Nodi. Utilizzando queste:

1) $ (z^(i+1)-2z^(i)+z^(i-1))/(p^2)-(z^(i+1)-z^(i-1))/(2p)*1/r^i-z^i/(r^i)^2=12/h^2*(z^(i+1)-z^(i-1))/(2p)*[ (u^(i+1)-u^(i-1))/(2p)+v*u^i/r^i+1/2*(z^i)^2 ]+qr^i/(2D) $

2) $ (u^(i+1)-2u^(i)+u^(i-1))/(p^2)-(u^(i+1)-u^(i-1))/(2p)*1/r^i-u^i/(r^i)^2=-(1-v)/(2*r^i)*(z^i)^2-z^i*(z^(i+1)-z^(i-1))/(2p) $

e le condizioni al contorno che sono:
$ u(a)=0 $
$ u'(b)=0 $
$ z(a)=0 $
$ z(b)=0 $

Ottengo il seguente sistema:

1) $ (z^(3)-2z^(2)+z^(1))/(p^2)-(z^(3)-z^(1))/(2p)*1/r^2-z^2/(r^2)^2=... $ Nodo 2 Equazione 1
2) $ (u^(3)-2u^(2)+u^(1))/(p^2)-(u^(3)-u^(1))/(2p)*1/r^2-u^2/(r^2)^2=... $ Nodo 2 Equazione 2
3) $ (z^(4)-2z^(3)+z^(2))/(p^2)-(z^(4)-z^(2))/(2p)*1/r^3-z^3/(r^3)^2=... $ Nodo 3 Equazione 1
4) $ (u^(4)-2u^(3)+u^(2))/(p^2)-(u^(4)-u^(2))/(2p)*1/r^3-u^3/(r^3)^2=... $ Nodo 3 Equazione 2
5) $ u^(1)=0 $ Condizione al contorno $ u(a)=0 $
6) $ 1/(2p)*[ 3*u^4-4*u^3+u^2 ] ]=0 $ Condizione al contorno $ u'(b)=0 $ (ho usato la formula per la differenza finita centrata al bordo)
7) $ z^1=0 $ Condizione al contorno $ z(a)=0 $
8) $ z^4=0 $ Condizione al contorno $ z(b)=0 $

Ci sono 8 equazioni e le 8 incognite sono $ z(1), z(2), z(3), z(4), u(1), u(2), u(3), u(4) $ . Il tutto forma un sistema non lineare di 8 equazioni in 8 incognite.
E' corretto?
Grazie mille veramente per l'aiuto
Zurtolo
Starting Member
Starting Member
 
Messaggio: 6 di 14
Iscritto il: 21/06/2019, 15:32

Prossimo

Torna a Analisi Numerica e Ricerca Operativa

Chi c’è in linea

Visitano il forum: Nessuno e 1 ospite