Euler implicite

hello, j'ai une ébauche de programme pour euler implicite mai il marche pas donc si quelqu'un pouvait joindre son programme ca serait sympa merci d'avance.
gab

Le prog avait déja été posté sur le forum ya qq temps. Le voila

Code:

clear
a=0; b=20; h=0.5; M=10;           
N=(b-a)/h;
x(1)=a;
y(1)=a;
g=inline('z-y-h*(1-(z/(2+sin(x))))');
j=inline('1+h*(1/(2+sin(x)))');
for k=2:N
    z(1)=y(k-1);
    for m=1:M
        z(m+1)=z(m)-(g(h,x(k-1),y(k-1),z(m))/j(h,x(k-1)));
    end
    y(k)=z(m+1);
    x(k)=x(k-1)+h;
end
plot(x,y);

merci simon

mais je comprend pas l'histoire du z qui prend la valeur du y et inversement dans le for du k dc si quelqu'un plus éclairé pouvait me répondre ca serait sympa merci d'avance. gab

y(n+1)=y(n) +y'(n+1)*h % c'est la déf d'euler implicite
=> y(n+1)-h*f(x(n+1),y(n+1))-y(n)=0=F(y(n+1))

Simon a décrété que y(n+1)=z et que y(n)=y. Puis il a remplacé dans F(y(n+1))=g
Par contre je cherche depuis un certain temps comment il est arrivé à j.

Peut-etre par Newton: x(N+1)=x(n)-f(x(n))/f '(x(n))

Je suis désolé ce programme n'est pas de moi, je l'ai simplement copié d'un autre topic dans le forum. Je peux pas te dire d'où vient le j. Je comprends pas très bien non plus.

j c'est la dérivée de g (par rapport a z),
et en effet, on utilise Newton pour cacluler z ( qui est en fait y(n) c-à-d l'inconnue dans l'algorithme d'Euler implicite )

Je ne comprend toujours pas cette partie-là...

Simon a écrit:

Code:


for k=2:N
    z(1)=y(k-1);
    for m=1:M
        z(m+1)=z(m)-(g(h,x(k-1),y(k-1),z(m))/j(h,x(k-1)));
    end
    y(k)=z(m+1);
    x(k)=x(k-1)+h;
end


à la base, on avait dit que z=y(n+1), et là avec les k jm'en sors plus.

et pi, pour z(m-1) là jsais d'où ça vient. mais (g(h,x(k-1),y(k-1),z(m))?? pourquoi x(k-1) ??

j'ai franchement dû louper qqch d'important là... o_O

merci:)

Ai envoyé un mail à Mr Tolley m'aider à comprendre Euler implicite, et je lui ai envoyé le programme qui circule sur le forum. il a répondu que le programme était foireux.

Il m'a donc envoyé une note, un programme et une fonction qui permettent de résoudre un cas particulier.

Je ne sais pas comment ajouter un .doc ici alors je vais devoir le mettre direct dans le post...

Notes données par Mr Tolley

1) Méthode de Newton

Elle sert à résoudre une équation à une inconnue du type f(x)=0.
L’inconnue est un réel (ou même un complexe).
Soit r le réel que l’on cherche. Si f(r)=0 mais f’(r)0 la racine de l’équation est simple et la méthode de Newton converge de manière quadratique. Si l’on prend comme valeur de départ du processus itératif, un x0 assez près de r on trouvera en très peu d’itérations (4 ou 5) une valeur très proche de celle de r (le nombre de chiffres précis double à chaque itération).

2) Méthode d’Euler explicite

Elle sert à résoudre un problème de Cauchy pour une équation ou un système d’équations différentielles ordinaires du premier ordre.

Dans le cas d’une seule équation le problème s’écrit :

y’(x) = f(x,y) et y(x0) = y0 a  x  b (évidemment : a = x0 et f(x,y) est donnée).

On cherche donc une fonction d’une variable réelle sur un intervalle donné [a,b] qui satisfasse à la fois à l’équation et à la condition initiale.

Sous certaines conditions qui concernent f(x,y) la solution de ce problème est unique.

L’algorithme d’Euler explicite consiste à approcher y(x) sur [a,b] par une suite de valeurs numériques à l’aide du processus itératif :

On pose x1 = x0 (= a) et y1 = y0
On choisit un nombre de points N pour partitionner l’intervalle [a,b] et on lui associe le pas d’intégration h = (b-a)/N.
On calcule les approximations de y aux différentes abscisses x(j) = a + j.h avec j = 1,2,…,N à l’aide du processus itératif :

x(1) = a ; y(1) = y0 ;
for j=1 : N
y(j+1) = y(j) + h*f(x(j),y(j)) ;
x(j+1) = a + j*h ;
end

Donc à chaque étape du processus on calcule un nombre, y(j+1) en fonction de quantités connues [y(j) et x(j) qui on été définis avant (pour j=1) ou calculés à l’étape précédente (pour j1 ; h défini avant et f(x(j),y(j)) que l’on calcule au moyen d’une fonction « inline » de deux variables que l’on doit aussi définir avant d’appliquer l’algorithme].

3) Méthode d’Euler implicite

Ici, les choses sont un peu plus subtiles. L’algorithme est différent et l’on devrait écrire :

x(1) = a ; y(1) = y0 ;
for j=1 :N
y(j+1) = y(j) + h*f(x(j+1),y(j+1)) ;
x(j+1) = a + j*h ;
end

Pour x(j+1) il n’y a pas de problème, il suffit de le définir avant de calculer y(j+1) [NB on peut même définir tous les x(j) avant la boucle].

Le hic, c’est qu’on cherche y(j+1) et que la formule qui fournit ce nombre contient le nombre lui-même !

On dit que y(j+1) est défini implicitement. On doit donc résoudre une équation à une inconnue que l’on peut écrire z = y(j) + h*f(x(j),z) où tout est connu sauf z qui n’est rien d’autre que y(j+1) .

Le truc consiste alors à résoudre z = y(j) + h*f(x(j),z) par la méthode de Newton en prenant comme z0 la valeur de y(j) calculée auparavant. Et ça marche bien si h est suffisamment petit.

En effet, si h est petit, y(j+1) ne sera pas trop différent de y(j) et si le problème de Cauchy n’admet qu’une solution (ce que l’on supposera), la méthode de Newton pour trouver y(j+1) convergera de manière quadratique, si bien qu’il suffira de faire un maximum de 5 itérations pour trouver y(j+1).

Code:

function E=eulim10_04(t0,tf,y0,N)
%
% Résout le problème de Cauchy
%
% y'(t) = -5/y² + t ;  y(0) = y0
%
% par la méthode d'Euler implicite à pas constant
%
% ENTREES
%
% t0 et tf sont les bornes de l'intervalle d'intégration
% y0 est la condition initiale : y(t0) = y0
% N est le nombre de pas d'intégration
%
% SORTIES
%
% E = [T' Y'] où T est le vecteurs des abscisses de calcul
% et Y celui des valeurs approchées de la fonction cherchée

h=(tf-t0)/N;
T=zeros(1,N+1); Y=zeros(1,N+1);
T=t0:h:tf;
Y(1)=y0;
% Application de la méthode de Newton pout trouver Y(j+1)
for j=1:N
    z(1)=Y(j);
    for i=1:10
        z(i+1)=z(i)-(z(i)^2*(z(i)-(Y(j)+h*(T(j)+h)))+5*h)/(3*z(i)^2-2*z(i)*(Y(j)+h*(T(j)+h)));
        if abs(z(i+1)-z(i)) < 1e-12
            v=z(i+1);
            break
        end
    end
    Y(j+1) = v;
end
E = [T' Y'];


% Résolution du problème de Cauchy
%
% y'(t) = -5/y² + t ;  y(0) = y0
%
% par la méthode d'Euler implicite (ou régressive).
%
% Appelle eulim10_04
%

H=[0.5 0.4 0.3 0.2 0.1 0.05];
t0=0; tf=15; y0=-1;
N=round((tf-t0)./H)+1;
%
figure(1); hold on; grid on;
axis([0 15 -3.5 0]); xlabel('t'); ylabel('y(t)')
titre=strcat('Problème de Cauchy      y''(t) = -5 / y²  +  t ;  y(0) = ',num2str(y0));
title(titre)
col=['k' 'r' 'b' 'm' 'k' 'r'];
for i=1:length(H)
E=eulim10_04(t0,tf,y0,N(i));
plot(E(:,1),E(:,2),col(i),'LineWidth',1)
pause
end
%
E=eulim10_04(t0,tf,y0,5000);
plot(E(:,1),E(:,2),'c','LineWidth',1.5)

arf arf arf... ce programme est d'une sombre obscurité ! il me semblait avoir plus ou moins compris euler implicite, ou du moind etre capable de résoudre un probleme avec mais la... je me remet en question!

- pq il fait T=zeros(1,N+1); Y=zeros(1,N+1);

- l'équation pour trouver z avec newton est étrange
z(i+1)=z(i)-(z(i)^2*(z(i)-(Y(j)+h*(T(j)+h)))+5*h)/(3*z(i)^2-2*z(i)*(Y(j)+h*(T(j)+h))); d'oû sort ce z(i)^2 et ce 5*h et ce 3*z(i)^2 ???
ca correspond en rien avec l'équation de départ, y a quelques éléments qu'on retrouve mais tellement peu, et ya surtout plin de truc bizares!

- est ce qu'il faut tjs modifer notre équation pour obtenir
g(z)=z-y-h*f(x(n+1),z)? ou bien c'est dla fumisterie ce truc?

bon j'espere que des lecteurs mieux avisés que moi auront limpidement compris ces nouvelles subtilités et qu'ils auront la gentillesse de les expliquer...
merci d'avance!

waw pq c'est écrit en tt petit, le monde est contre moi ajd hui!!
ooo désespoir! jvoudrais etre de retour a l'age de pierre et devoir suer pour me nourir plutot que de me pisser dessus devant cet ordi de malheur!!

Résolution du problème de Cauchy par Méthode de Euler implicite:
Euler implicite : y(n+1)=y(n) +y'(n+1)*h
Newton: y(x(n+1))=y(x(n))-g(x(n))/g’(x(n))
y'(t)=f(x,y)= x-5/(y.^2)
=> y(n+1)-h*f(x(n+1),y(n+1))-y(n)=0=g(z)
z=y(n+1); x(n+1)=x(n)+h; f(x(n+1),y(n+1))=x(n)+h-5/(z(i).^2);

il faut remplacer tout ça dans l'eq de Newton:

Newton: z(i+1)=z(i)- g(x(i+1),y(i+1))/g'(x(n+1),y(n+1))


z(i+1)=z(i)-(z(i)^2*(z(i)-(Y(n)+h*(T(n)+h)))+5*h)/(3*z(i)^2-2*z(i)*(Y(n)+h*(T(n)+h)));

et T=zeros(1,N+1);Y=zeros(1,N+1); sont 2 matrices à 1 ligne et N+1 colonnes ! (ce sont en fait 2 vecteurs représentés respectivement sur les axes X et Y)

Ben en gros c'est comme on fait nous sauf que lui à tout décomposer à l'intérieur de la formule... non ?

enfait, la version qui est sur le forum, je ne la comprend pas, donc je ne sais pas répondre à ta question.
Je dois t'avouer que je préfère la version de Mr Tolley, car il avait sérieusement insisté sur le fait que le prog était foireux et qu'il a même pas essayé de le comprendre...

Ca m'a l'air bien compliqué tout ces bazars... voila une version que je suis sur qui fonctionne... j'avoue j'ai rajouté un truc pour faire varier le pas et chager la couleur de la courbe selon le pas, mais l'intérieur du programme est l'essentiel (ce qui est contenu dans le premier 'for')

clear
clc

col=['r' 'b' 'g' 'c' 'm' 'y'];
h=[0.5 0.2 0.1 0.05 0.01 0.005];
g = inline ('z-y-h*(1-(z/2+sin(x)))');
gp = inline ('1+h*(1/2+sin(x))');
for k=1:6
a=0;
b=20;
N=round((b-a)/h(k));
x(1)=a;
y(1)=0;
for n=2:N
x(n)=h(k)+x(n-1);
z(1)=y(n-1);
for m=2:10
z(m)=z(m-1)-(g(h(k),x(n-1),y(n-1),z(m-1))/gp(h(k),x(n-1)));
end
y(n)=z(10);
end
plot (x,y,col(k))
hold on
grid on
axis equal
pause
end

c'est pas exactement la meme chose que lautre qui est sur le forum ca?

C quoi cette histoire de round??
dans N= round......

c'est une fonction matlab qui prend l'entier le plus proche. Comme N c'est le nombre d'itérations, ça serait gênant d'avoir 17.2 itérations a faire . Bon j'avoue avec les h qui sont la je pense pas qui ait de chance de tomber sur autre chose qu'un entier, donc le round est totalement inutile ...