A.A. 2019/20
docente: Sandra Carillo
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;
ivp := ode({f''(x) + 2*0.1*f'(x) + f(x)=0, f(0) = 1, f'(0) = 0}, f(x))
la cui soluzione
esatta,indicata
con (3), e` data da:
Fes(x):=solve(ivp);
Fex(x):=simplify(%)
che possiamo anche scrivere come:
F(x):=Fex(x);
indicata con (5), il cui
grafico, e`
plotfunc2d(exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)),x=0..PI)
nell'intervallo 0..PI, e
plotfunc2d(exp(-.1*x)*(cos(0.9949874371*x)+0.1005*sin(0.9949874371*x)),x=
0..7*PI)
f(x,`ϵ`):=f_0(x) +`ϵ`*f_1(x) +`ϵ`^2 *f_2(x) +`ϵ`^3* f_3(x) +`ϵ`^4 * f_4(x);
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);
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);
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,
f_0(0)=1;
e (11) agli ordini
successivi, cioe` per k=1,2, ...
f_k(0)=0;
ed il valore iniziale della
derivata prima della funzione nel punto x=0
=b
f_0'(0)=0;
e le condizioni (13), agli ordini
successivi, cioe` per k=1,2, ...
f_k'(0)=0;
der2+ 2*`ϵ`*der1+f(x,`ϵ`)=0;
ode2:=%:
Passo 2
Soluzione del problema all'ordine
zero: consideriamo, cioe`, epsilon=0 nella (14):
subs(ode2,`ϵ`=0);
delete f0:
ivp0 := ode({f_0''(x) + f_0(x)= 0, f_0(0) = 1, f_0'(0) = 0}, f_0(x))
f0(x):=solve(ivp0);
f0_primo(x):=expand(diff(cos(x),x));
ode2;
expand(ode2,`ϵ`);
delete f1;f0(x);
ivp1 := ode({f1''(x) + f1(x)= 2*(sin(x)), f1(0) = 0, f1'(0) = 0}, f1(x))
f1(x):=solve(ivp1);
In generale, otteniamo:
f(x)=f0(x)+`ϵ`*f1(x);
che, nel caso a) diventa
f(x)=f0(x)+0.1*f1(x);
plotfunc2d(cos(x)+0.1*(sin(x)-x*cos(x)), x= 0..PI)
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)
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)
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)
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/`ϵ`);
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 sistemache 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;
c) (attenzione 0.5 non e` << 1 !) `ϵ`=.5;
b)
ivp := ode({f''(x) + 2*0.01*f'(x) +f(x)=0, f(0) = 1, f'(0) = 0}, f(x))
la cui soluzione esatta,indicata con (3), e` data da:
Fes(x):=solve(ivp);
nella quale mettiamo in evidenza il termine esponenziale, ottenendo la
seguente, indicata con (4):
Fex(x):=simplify(%)
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)
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)
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)
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)
ivpk := ode({fk''(x) + fk(x)=-2*(diff(f(k-1)(x), x)), fk(0) = 0, fk'(0) = 0}, fk(x));
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))
Fes(x):=solve(ivp);
Fex(x):=simplify(%)
il cui grafico e`
plotfunc2d(exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)),x= 0..PI)
plotfunc2d(exp(-.01*x)*(cos(0.9999499987*x)+0.1000050004*sin(0.9999499987*x)),x= 0..7*PI)
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));
plotfunc2d(cos(x)+.01*(sin(x)-x*cos(x)), x= 0..PI)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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
Ultimo caso c) `ϵ`=.5;
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));
float(sqrt(3)/2);
f:=expand(solve(ivp));
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)
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