-
Notifications
You must be signed in to change notification settings - Fork 0
/
eulerexpl.m
57 lines (48 loc) · 1.39 KB
/
eulerexpl.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
function [t,y] = eulerexpl(f,t0,y0,h,nbpas)
%
% MŽethode d'Euler explicite pour les syste?mes d'eŽquations differentielles
% de la forme y'(t) = f(t,y(t)) {SD}
%
% PrŽealable:
% Vous devez creer un fichier .m contenant la fonction f(t,y).
% Voir par exemple le fichier eqndiff.m
%
% Exemple d'appel :
% [t,y] = eulerexpl('eqndiff',0.0,[1 1],1.0e-1,10)
%
% Arguments en entrŽe:
% 1) f: Le nom entre apostrophes (' ') du fichier .m contenant la fonction
% f(t,y).
% 2) t0: le temps initial
% 3) y0: la (les) condition(s) initiale(s), sous la forme d'un vecteur ligne :
% par exemple, [1 2]. La dimension de y0 donne la dimension du syst?eme
% 4) h: le pas de temps
% 5) nbpas: le nombre maximal de pas de temps
%
% Arguments en sortie :
% 1) t est un vecteur contenant les pas de temps, ti
% 2) y est une matrice contenant la solution obtenue :
% la colonne i de y correspond a la solution de l'eŽquation i
% on fait quelques calculs prŽeliminaires
nbeq = size(y0,2);
y0 = y0';
k1 = zeros(nbeq,1);
y = zeros(1,nbeq);
n = 1;
t(1) = t0;
y(1,:) = y0';
warning off
% la boucle transitoire
while(n <= nbpas),
% on fait une couple de vŽrification
k1 = feval(f,t(n),y0);
if (any(~isfinite(k1))),
disp('la fonction n''est pas definie en certains points');
return;
end
% la mŽthode d'Euler est ici!
y0 = y0 + h * k1;
y(n+1,:) = y0';
t(n+1) = t(n) + h;
n = n+1;
end