A.A. 2019/20
docente: Sandra Carillo

Esercitazione 1
Metodo Perturbativo Diretto

con l'ausilio del programma di calcolo simbolico
ToolBox in MatLab

MuPAD http://www.mupad.de/


Problema dell'oscillatore debolmente smorzato


 

In questa esercitazione consideriamo il problema
a
valori iniziali seguente,indicato con (1):




dove l'incognita e` la funzione

ivp := ode({f''(x) +2*`ϵ`*f'(x) + f(x)= 0, f(0) = 1, f'(0) = 0}, f(x))


Vogliamo, cioe`, studiare il moto di un punto materiale di massa m=1, soggetto ad una forza elastica, con costante elastica K=1,
vincolato bilateralmente ad una linea orizzontale in presenza di un attrito
piccolo.

Consideriamo i casi


a)
   `ϵ`=.1;



b)   `ϵ`=.01;




c) (attenzione 0.5 non e` piccolo !)  `ϵ`=.5;



Caso a) sostituendo il valore di epsilon nella (1), otteniamo il seguente ivp,indicato con (2):

ivp := ode({f''(x) + 2*0.1*f'(x) + f(x)=0, f(0) = 1, f'(0) = 0}, f(x))

math

la cui soluzione esatta,indicata con (3), e` data da:

Fes(x):=solve(ivp);

math

nella quale mettiamo in evidenza il termine esponenziale, ottenendo la seguente, indicata con (4):

Fex(x):=simplify(%)

math

che possiamo anche scrivere come:

F(x):=Fex(x);

math

indicata con (5), il cui grafico, e`

plotfunc2d(exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)),x=0..PI)

MuPAD graphics

nell'intervallo 0..PI, e

 

plotfunc2d(exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)),x=

0..7*PI
)

MuPAD graphics

nell'intervallo 0..7*PI. Adesso calcoliamo la soluzione mediante lo sviluppo perturbativo diretto,
cioe` sostutiamo (6)

 

 

f(x,`ϵ`):=f_0(x) +`ϵ`*f_1(x) +`ϵ`^2 *f_2(x) +`ϵ`^3* f_3(x) +`ϵ`^4 * f_4(x);

math

 

le cui derivate forniscono, rispettivamente, (7):


 

 

der1:= diff(f_0(x),x) +`ϵ`*diff(f_1(x),x) +`ϵ`^2 *diff(f_2(x),x) +`ϵ`^3 * diff(f_3(x),x) +`ϵ`^4 * diff(f_4(x),x);

math

 

e (8):


 

 

der2:= diff(diff(f_0(x),x),x)+`ϵ`*diff(diff(f_1(x),x),x) +`ϵ`^2 *diff(diff(f_2(x),x),x) +
`ϵ`^3 *
diff(diff(f_3(x),x),x) +`ϵ`^4 *  diff(diff(f_4(x),x),x);

math

 

ed imponiamo le condizioni iniziali, indicate con (9)


=a

 

 

nel caso che stiamo cosiderando a=1 e` il valore iniziale della funzione nel punto x=0, quindi,
all'ordine zero in `epsilon' imponiamo (10):


f_0(0)=1;

math

e (11) agli ordini successivi, cioe` per k=1,2, ...

f_k(0)=0;

 

math

ed il valore iniziale della derivata prima della funzione nel punto x=0

=b

nel caso che stiamo cosiderando b=0 e` il valore iniziale della derivata prima della funzione incognita nel punto x=0,
quindi, all'ordine zero in epsilon imponiamo (12):

f_0'(0)=0;

math

e le condizioni (13), agli ordini successivi,  cioe` per k=1,2, ...

f_k'(0)=0;



L'equazione differenziale dell'oscillatore smorzato, nella quale siano state sostituite le (6), (7) e (8), assume la forma (14):

 

der2+ 2*`ϵ`*der1+f(x,`ϵ`)=0;


 

ode2:=%:

                                     Passo 2

Soluzione del problema all'ordine zero: consideriamo, cioe`, epsilon=0 nella (14):

subs(ode2,`ϵ`=0);


che, considerando le condizioni inziali (12) e (13), fornisce il seguente problema di Cauchy, (15)

delete f0:

ivp0 := ode({f_0''(x) + f_0(x)= 0, f_0(0) = 1, f_0'(0) = 0}, f_0(x))

math

f0(x):=solve(ivp0);

math

 

f0_primo(x):=expand(diff(cos(x),x));

math

ode2;

math

expand(ode2,`ϵ`);

 

math

 

math

i termini del primo ordine in epsilon, con le condizioni iniziali omogenee per f1, cioe` f1(0)=f1'(0)=0, forniscono il seguente problema di Cauchy, (16)

delete f1;f0(x);


 

ivp1 := ode({f1''(x) +  f1(x)= 2*(sin(x)), f1(0) = 0, f1'(0) = 0}, f1(x))

math

f1(x):=solve(ivp1);

math

In generale, otteniamo:

f(x)=f0(x)+`ϵ`*f1(x);

math

 

che, nel caso a) diventa

f(x)=f0(x)+0.1*f1(x);

math

                             Passo 3
 

Stampa dei grafici della funzione f(x) approssimata al primo ordine, in vari intervalli:

plotfunc2d(cos(x)+0.1*(sin(x)-x*cos(x)), x= 0..PI)

MuPAD graphics


ed il confronto, nello stesso intervallo, tra tale soluzione, in blu, e la soluzione esatta, in rosso, cioe`:

plotfunc2d(cos(x)+0.1*(sin(x)-x*cos(x)),exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)),

x= 0..PI
)

MuPAD graphics


Come si vede, l'accordo tra le due soluzioni appare tanto migliore quanto piu` x e` prossimo a zero.
Scegliamo ora l'intervallo [0,3 PI]

plotfunc2d(cos(x)+0.1*(sin(x)-x*cos(x)),exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)),

x= 0..3*PI
)

MuPAD graphics

e poi l'intervallo [0,9 PI]

plotfunc2d(cos(x)+0.1*(sin(x)-x*cos(x)),exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)),x=0..9*PI)

 

MuPAD graphics

Si osserva che, per x tale che

O(`ϵ`*x)=1; cioe`,


 



cioe` x e` di ordine 1/epsilon e, quindi, il termine f_1(x) e` dello stesso ordine del termine

f_0(0) e, quindi, non rappresenta una `piccola' correzione, corrispondentemente, osservo la crescita lineare della f_1(x) che non verifica la condizione di limitatezza



sotto la quale lo sviluppo perturbativo e` valido


x=O(1/`ϵ`);

math

la soluzione approssimata ha un comportamento completamente diverso da quella esatta.

N.B.la soluzione approssimata ha ampiezza che cresce linearmente con x: inaccettabile poiche` questo corrisponde ad energia crescente contrariamente al fatto che il sistema

che si sta considerando e` dissipativo (anche se debolmente).

 

Infine, confrontiamo la soluzione approssimata con quella esatta in corrispondenza

ai due valori di epsilon nei casi


b) `ϵ`=.01;

math

c) (attenzione 0.5 non e` << 1 !)  `&epsiv;`=.5;

math

b)

 

ivp := ode({f''(x) + 2*0.01*f'(x) +f(x)=0, f(0) = 1, f'(0) = 0}, f(x))

math

la cui soluzione esatta,indicata con (3), e` data da:

Fes(x):=solve(ivp);

math

nella quale mettiamo in evidenza il termine esponenziale, ottenendo la

seguente, indicata con (4):

Fex(x):=simplify(%)

math

 

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)),x= 0..3*PI)

MuPAD graphics

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)), x= 0..9*PI)

MuPAD graphics

plotfunc2d(cos(x)-0.1*x*cos(x),exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)), x= 0..15*PI)

MuPAD graphics

plotfunc2d(cos(x)-0.1*x*cos(x),exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)), x= 0..9*PI)

 

MuPAD graphics

Naturalmente, un modo per migliorare l'approssimazione e` quello di calcolare termini successivi, cioe` non fermarsi all'ordine k=1, ma risolvere altri problemi: tutti del tipo:

 

ivpk := ode({fk''(x) + fk(x)=-2*(diff(f(k-1)(x), x)), fk(0) = 0, fk'(0) = 0}, fk(x));

math

cioe` lineari (stesso operatore differenziale), con condizioni iniziali omogenee, non omogenei, nei quali il termine noto (forzante) e` legato alla soluzione del problema differenziale al passo precedente.

 

D'altra parte vediamo che l'accordo tra la soluzione esatta e quella ottenuta mediante il metodo perturbativo diretto (al primo ordine) quanto piu` x e` prossimo a zero. Adesso vediamo qual'e` l'andamento quando scegliamo

 

b) epsilon=.01

 

ivp := ode({f''(x) + 2*0.01*f'(x)+f(x)=0, f(0) = 1, f'(0) = 0}, f(x))

math

 

Fes(x):=solve(ivp);

math

 

Fex(x):=simplify(%)

math

il cui grafico e`

plotfunc2d(exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)),x= 0..PI)

MuPAD graphics

 

plotfunc2d(exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)),x= 0..7*PI)

MuPAD graphics

nell'intervallo 0..7*PI. Adesso confrontiamo tale soluzione con quella ottenuta mediante lo sviluppo perturbativo diretto.


Osservazione: non abbiamo bisogno di calcolare nuovamente la soluzione perche` i problemi ad ogni ordine sono indipendenti dal valore di epsilon.

Quindi, la soluzione approssimata al primo ordine e`:

f(x,0.1)=cos(x)+.01*(sin(x)-x*cos(x));

math

 

plotfunc2d(cos(x)+.01*(sin(x)-x*cos(x)), x= 0..PI)

MuPAD graphics

Ora confrontiamo l'andamento della soluzione esatta (in rosso) con quella approssimata al primo ordine (in blu) nell'intervallo [0,PI] nel grafico seguente

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)),x= 0..PI)

MuPAD graphics

nell'intervallo [0,3 PI]  nel grafico seguente

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)),x= 0..3*PI)

MuPAD graphics

nell'intervallo [0,10 PI] nel grafico seguente

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)),x= 0..10*PI)

MuPAD graphics

nell'intervallo [0,20 PI] nel grafico seguente

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)), x= 0..20*PI)

MuPAD graphics

nell'intervallo [0,40 PI] nel grafico seguente

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)), x= 0..40*PI)

 

MuPAD graphics

in quest'ultimo grafico si raggiungono valori di x per i quali la soluzione al primo ordine (blu) non e` accettabile.

Confronto tra le soluzioni con i due diversi valori di epsilon

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),cos(x)+0.1*(sin(x)-x*cos(x)),exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)), exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)),x=0..2*PI)

 

MuPAD graphics

nell'intervallo [0,2 PI] le curve rossa e verde rappresentano,rispettivamente, le soluzioni approssimata ed esatta per epsilon=0.1 mentre le curve blu e senape,rispettivamente, le soluzioni approssimata ed esatta per epsilon=0.01

 

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),cos(x)+0.1*(sin(x)-x*cos(x)),exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)), exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)), x=0..PI)

 

MuPAD graphics

nell'intervallo [0, PI] le curve rossa e verde rappresentano<>,rispettivamente, le soluzioni approssimata ed esatta per epsilon=0.1 mentre le curve blu e senape,rispettivamente, le<> soluzioni approssimata ed esatta per epsilon=0.01

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),cos(x)+0.1*(sin(x)-x*cos(x)),exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)), exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)), x=0..3*PI)

 

 

MuPAD graphics

nell'intervallo [0,3 PI] le curve rossa e verde rappresentano,rispettivamente, le soluzioni approssimata ed esatta per epsilon=0.1 mentre le curve blu e senape, rispettivamente, le soluzioni approssimata ed esatta per epsilon=0.01

plotfunc2d(cos(x)+0.01*(sin(x)-x*cos(x)),cos(x)+0.1*(sin(x)-x*cos(x)),exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)), exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)), x=0..6*PI)

 

MuPAD graphics

nell'intervallo [0,6 PI] le curve rossa e verde rappresentano,rispettivamente, le soluzioni approssimata ed esatta per epsilon=0.1 mentre le curve blu e senape,rispettivamente, le soluzioni approssimata ed esatta per epsilon=0.01.

Si nota, dal confronto tra il grafico della soluzione approssimata (blu) e quello della corrispondente soluzione esatta (senape), che, quando epsilon=0.01, la differenza tra le due soluzioni e` percentualmente limitato. Al contrario, dal confronto tra il grafico della soluzione approssimata (rossa) e quello della corrispondente soluzione esatta (verde), che, quando epsilon=0.1, la differenza tra le due soluzioni e` crescente al variare di x ed e` di ordine 1 per
 image

Ultimo caso c)  `&epsiv;`=.5;

math
Innanzitutto, devo calcolare la soluzione esatta

 

delete f:

 

ivp := ode({f''(x) + 2*0.5*f'(x) + f(x)=0, f(0) = 1, f'(0) = 0}, f(x));

math

 

float(sqrt(3)/2);

math

f:=expand(solve(ivp));

math

Invece, per ottenere quella approssimata, devo solo sostituire ad epsilon il valore 0.5.
Ottengo

 

plotfunc2d(cos(x)+0.5*(sin(x)-x*cos(x)),exp(-.5*x)*(cos((sqrt(3)/2)*x)+(sqrt(3)/3)*sin((sqrt(3)/2)*x)), x= 0..PI)

MuPAD graphics

Come previsto, dal confronto tra il grafico della soluzione approssimata (blu) e quello della corrispondente soluzione esatta (rossa), che, quando epsilon=0.5, la differenza tra le due soluzioni e` percentualmente elavata gia` per x dell'ordine di PI/2. D'altra parte, come gia` osservato, .5 non e` << 1