Lezione del 9 marzo 2007¶
Metodi multistep e instabilità
Funzione di test:
File fun_exp.m
:open:
1%
2% funzione standard per testare la stabilità
3%
4function yp=fun_exp(x,y)
5 yp=-10*y;
6end
Metodo multistep instabile:
File instabile.m
:open:
1%
2% Risolve la ODE
3% y' = -a*y
4% y(0) = 1
5%
6% con il metodo
7% y(1) = y(0) + h f(x(0),y(0))
8% y(k+1) = y(k-1) + 2 h f(x(k),y(k))
9%
10function [x,y]=instabile(fun,xf,N)
11 h = xf / N ;
12 x(1) = 0 ;
13 y(1,:) = 1 ;
14 %% passo Eulero Esplicito
15 x(2) = x(1) + h ;
16 f1 = feval(fun,x(1),y(1,:));
17 y(2,:) = y(1,:) + h * f1 ;
18 %% ciclo con il metodo multistep
19 for k=2:N
20 fk = feval(fun,x(k),y(k,:));
21 x(k+1) = x(k)+h ;
22 y(k+1,:) = y(k-1,:)+2*h*fk ;
23 end
24end
Programma driver per verificare il codice:
File test_instabile.m
:open:
1N = 100 ;
2[x,y] = instabile (@fun_exp,1,N) ;
3hold off ;
4plot(x,y,'.-b');
5hold on ;
6x = [0:0.01:1];
7y = exp(-10.*x) ;
8plot(x,y,'-r','LineWidth',2);
9axis([0 1 -1 2]);
10legend('Instabile','Esatta');
Esempio di output:
Metodo multistep instabile con condizioni iniziali esatte:
File instabile1.m
:open:
1%
2% Risolve la ODE
3% y' = -a*y
4% y(0) = 1
5%
6% con il metodo
7% y(1) = 1/(sqrt(1+alpha^2*h^2) + h*alpha)
8% y(k+1) = y(k-1) + 2 h f(x(k),y(k))
9%
10function [x,y]=instabile1(xf,N)
11 h = xf / N ;
12 x(1) = 0 ;
13 y(1,:) = 1 ;
14 alpha = 10 ;
15 %% passo Eulero Esplicito
16 x(2) = x(1) + h ;
17 y(2,:) = 1/(sqrt(1+alpha^2*h^2) + h*alpha) ;
18 %% ciclo con il metodo multistep
19 for k=2:N
20 fk = - alpha * y(k,:) ;
21 x(k+1) = x(k)+h ;
22 y(k+1,:) = y(k-1,:)+2*h*fk ;
23 end
24end
Programma driver per verificare il codice:
File test_instabile.m
:open:
1N = 100 ;
2[x,y] = instabile (@fun_exp,1,N) ;
3hold off ;
4plot(x,y,'.-b');
5hold on ;
6x = [0:0.01:1];
7y = exp(-10.*x) ;
8plot(x,y,'-r','LineWidth',2);
9axis([0 1 -1 2]);
10legend('Instabile','Esatta');
Esempio di output:
