Risolvere numericamente Sistema di equazioni differenziali con valori ai bordi
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 per (https://www.matematicamente.it/forum/vi ... p?t=189081)
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
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);
endCodice per (https://www.matematicamente.it/forum/vi ... p?t=189081)
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
Risposte
Ciao, un paio di precisazioni innanzitutto
Indipendente?
\(b < a\)?
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.
"Zurtolo":
* r è la variabile dipendente
Indipendente?
"Zurtolo":
L'integrazione deve avvenire sull'intervallo [a b], dove:
*a=12
*b=0
\(b < a\)?
"Zurtolo":
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.
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
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
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
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
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":
$ 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":
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":
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.
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
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
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.
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.
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
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":
*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":
**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":
***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.
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
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
Mi sembra corretto. Al posto dell'equazione 6 puoi più semplicemente imporre che \(u^3 = u^4\).
Grazie mille al rientro dalle ferie proverò a implementare il tutto e a vedere se funziona. Grazie ancora per l'aiuto, un saluto a tutti
Ciao e scusa il disturbo. Volevo chiederti se potevi condividere lo script attraverso il quale hai risolto il tuo sistema di equazioni differenziali con la funzione bvp4c perchè cercavo di risolvere un problema simile con tale funzione senza avere un esito positivo.
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo