Séance 11 - Corrigé

%seance 10 --> révision exercice 1
clc
clear
%point fixe

f=inline('-x+2-3*cos(5*x)')
fp=inline('-1+15*sin(5*x)')
g1=inline('2-3*cos(5*x)')
g2=inline('x-(-x+2-3*cos(5*x))/(-1+15*sin(5*1))')
x(1)=1;
u=0:0.01:10;
y=f(u);
plot(u,y); hold on; zoom on; grid on;
x(1)=1;
for n=1:15
x(n+1)=g2(x(n));
end
x'
pause
%newton
x(1)=1;
for n=1:5
x(n+1)=x(n)-f(x(n))/fp(x(n));
end
format long;
x'
pause
%dichotomie

a=1; b=1.2; sa=sign(f(a)); sb=sign(f(b));
for n=1:40
if sb==sa
break
disp('pu')% la valuer réelle est hors de l'intervalle!
else w=(a+b)/2;
sw=sign(f(w));
if sw==sa
a=w;
else
b=w;
end
end
xD(n)=w;
end
xD'

merci!
pour l exo 2 est ce que qqun sait quelle méthode utiliser ds quelles conditions?qqun aurait également un bon programme ac la méthode de heun?

%seance 11 exercice 2
%euler explicite
clc
clear
f=inline ('-5/y^2+sin(x)/5','x','y')%attention dans inline, il faut absolument donner les parametres dela foncctions s'il sont autres que x!!
y(1)=-1;
x(1)=0;
b=15; a=0;
h=0.05;
N=(b-a)/h;
for n=1:N
y(n+1)=y(n)+h*f(x(n),y(n));
x(n+1)=x(n)+h;
end
plot (x,y,'r'); grid on; zoom on; hold on;
% dans se prog,on cherche a approximer une fonction a partir de sa dérivée
%on obtient donc comme réponse un graph de fonction f(x)
pause

%euler implicite
g=inline('z-y-h*(-5/z^2+sin(x)/5)')
gp=inline('1+(10*h)/z^3')
for k=2:N%attention mettre 2 car sinon quand on fait y(n-1) si on avait pas 2 ca ferai 0 et matlab il aime pasca!!!
z(1)=y(k-1);
x(k)=x(k-1)+h;
for i=1:5% je fais newton pour trouver les y(n) qui sont inconnus!!
z(i+1)=z(i)-g(h,x(k-1),y(k-1),z(i))/gp(h,z(i));%je regarde dans mon inline tte les variables ou constantes qui interviennes, et il faut les définir, pour que quand j'appel le inline, il puisse savoir ce qu'il faut sur ces variables
end %la premiere boucle fais varier les x et y, donc dans la deuxieme boucle, x et y varient avec l'exposant du premier for.
y(k)=z(i)
end
plot(x,y,'b'); hold on; grid on;
pause
%heun

for n=2:N % attention a n=2!! sinon x(n-1) donnera x(0) et il est pas contant!
x(n)=x(n-1)+h;
k1=f(x(n-1),y(n-1));
k2=f(x(n),y(n-1)+h*k1);
y(n)=y(n-1)+h/2*(k1+k2);
end
plot (x,y'g'), hold on; grid on;

tout bon nemo merci boco

Salut les les gars,
Voici un petit exmple pour lexercice 3, je ne certifie pas lexactitude de cette chose, mais ça fonctionne...alors...alors voilà

%problème au conditions limites
%discrétisation
clc
clear
xdebut=0;xfin=pi;
N=100;
h=(xfin-xdebut)/N;
x=xdebut:h:xfin;
ydebut=1;
yfin=-1;

%écriture des matrice A et B
for n=1:length(x) %ici on ne fait pas length -1 car il nous faut 101 valeur
%tout comme le vecteur x
Ap(n,n)=-2+sin(x(n));
B(n,1)=(x(n)-1)*h^2;
end
A= Ap+diag(ones(N,1),1) + diag(ones(N,1),-1);

%il faut penser à utiliser les conditions initiales, lorsqu'on calcul B(1,1)
%on na pas tenu compte de y0 idem pour B(length(x)-1,1) on a pas tenu compte
%de yfinal
B(1,1)=B(1,1)-ydebut;
B(length(x),1)=B(length(x),1)-yfin;

Y=A\B
plot(x,Y')
grid on

Nemo, il ne faut tout de meme pas abuser de matlab... je veux dire que quand on lit tes messages ça fait peur!
on dirait presque que tu t es transformé en programme et quand on te pose une question tu répond en matlab!

va t aérer!
lol

hello, pour la méthode de dichotomie je sais que suivant la valeur du test f(a).f(w) selon qu'il est positif ou négatif le changement de signe aura lieu dans le premier intervalle ou dans le second ce qui situe le zéro dans le premier ou le second intervalle puisque le zéro correspond au changement de signe et si cette valeur vaut 0 on a trouvé notre zéro.
bon c bien beau tout ca mais ca me di pas grand chose pour matlab... j'ai bien remarqué la fct signe qui nous donne le signe de la fct dans le premier intervalle et dans le second donc selon a et selon B mais alors le reste...en fait je ne comprend pas à quoi correspond les égalités... merci d'avance pour votre réponse.

c ok pour la dichotomie il y a 5 lignes d'explication dans le résumé de verlant-chenet et ca a suffit...

Marie a écrit:Salut les les gars,
Voici un petit exmple pour lexercice 3, je ne certifie pas lexactitude de cette chose, mais ça fonctionne...alors...alors voilà

%problème au conditions limites
%discrétisation
clc
clear
xdebut=0;xfin=pi;
N=100;
h=(xfin-xdebut)/N;
x=xdebut:h:xfin;
ydebut=1;
yfin=-1;

%écriture des matrice A et B
for n=1:length(x) %ici on ne fait pas length -1 car il nous faut 101 valeur
%tout comme le vecteur x
Ap(n,n)=-2+sin(x(n));
B(n,1)=(x(n)-1)*h^2;
end
A= Ap+diag(ones(N,1),1) + diag(ones(N,1),-1);

%il faut penser à utiliser les conditions initiales, lorsqu'on calcul B(1,1)
%on na pas tenu compte de y0 idem pour B(length(x)-1,1) on a pas tenu compte
%de yfinal
B(1,1)=B(1,1)-ydebut;
B(length(x),1)=B(length(x),1)-yfin;

Y=A\B
plot(x,Y')
grid on

Hello, je crois que j'ai trouvé une petite faute dans ce programme, pour moi la matrice A devrait être remplie comme ceci :
Ap(n,n)=-2+(h^2)*sin(x(n));
Chez toi, il manque le facteur h^2 devant le sinus. Si tu remplaces y'' dans l'équation de départ pour trouver A et B, faut multiplier chaque membre de l'équation par h^2 pour le faire disparaitre du dénominateur de y''.

De plus, si on recherche sur ton graphe, les conditions limites, on les retrouve pas

Alors que si on rajoute le facteur manquant, ça marche.

Bon je met ma version de ce programme, mais j'ai deux petites questions :
1) ça sert a quelque chose de faire A=eye(N); ?
2) la méthode de tir au but, vue à la séance 9, c'est moi ou on a aucune théorie dessus ? Risque-t-on vraiment de l'avoir à l'exam?


Mon programme :

Code:

clc
clear

xfin=pi;
xdebut=0;
y0=1;
ypi=-1;

N=400;

h=(xfin-xdebut)/(N+1);

%A=eye(N); % je vois pas l'interet de ce truc... quelqu'un a-t-il une idée?

for n=1:N
    x(n)=h*n;
    A(n,n)=(-2+(h^2)*sin(x(n)));
    B(n)=(h^2)*(x(n)-1);
end

B(1)=B(1)-y0;
B(N)=B(N)-ypi;

A=A+diag(ones(N-1,1),1)+diag(ones(N-1,1),-1);

y=A\B';

plot(x,y,'g')
grid on;

pour le A=eye(N), ca sert à créer une matrice identité de taille N, donc avec des 1 sur la diagonal et des zéro aux autres endroits. du coup t'a plus ka multiplier cette matrice par la valeur que tu veux obtenir pour ta diagonale, ici en l'occurence -2+h^2*sin(x(n)).

et pour la théorie, y a rien a ce qu'il me semble dans ce ke le prof a donner et ki est sur le forum mais a la séance 9 l'assistante a donner des feuilles de théories. je pense ke la méthode des matrice est meilleure car elle fonctionne tjs et contrairement a la méthode du tir au but il faut pas tâtonner pour trouver des valeurs. dans le tir au but, il faut chipotter pour trouver C et c'est casse couille. bref personnellement j'vé plutot me concentré sur stafaire des matrice!!

ptite question pour cette méthode des matrice (qui doit bien avoir un nom un peu plus précis que ca j'imagine...) donc jvoudrais savoir comment on discrédirte un Y', au tp il donnait pour y''=y(n+1)-2y(n)+y(n-1)/h^2 ok c'est nikel pour cet exo mais si jamais on a du y'!!? qu'est ce qui ce pass?
si quelqu'un pouvait m'renseigner ce serai bien aimable de ca part!
merci

Voilà ce que j'ai pour l'exo 2 avec euler implicite
C'est pas la même chose que Nemo mais mon graph à l'air moins bizar que le sien... Donc si une troisième personne à qqch, ce serait bien de voir quel programme est bon

%Exo 2 par euler implicite
clear;
clc;
echo on
%Je vai choisir la méthode d'Euler implicite car parmi les méthodes qu'on a
%vues, c'est celle qui ne présente aucun risque d'instabilité
echo off
f=inline('-5/(y.^2)+sin(x/5)');
fprime=inline('10/(y.^3)');
a=0;
b=15;
N=200; %Je prend un n délibérément grand car au plus N est grand, au plus la précision de la réponse sera grande
h=(b-a)/N;
x(1)=a;
y(1)=-1;%Condition à l'origine
for k=1:N
x(k+1)=x(k)+h; %Création du vecteur x
%Normalement on aurait ensuite: y(k+1)=y(k)+h*f(x(k+1),y(k+1)); mais on ne peut
%pas appliquer cette formule car y(k+1) est inconnu.
%On a donc l'équation z-y(k)-h*f(x(k+1),z)=0 à résoudre en z, ce qui
%nous donnera z donc y(k+1). Je vai résoudre cette équation par la
%méthode de Newton.
z=y(k);
for n=1:10
%z(n+1)=z(n)-g(z(n))/gprime(z(n));
%g(z)=z-y(k)-h*f(x(k+1),z);gprime(z)=1-h*fprime(x(k+1),z)
z(n+1)=z(n)-(z(n)-y(k)-h*f(x(k+1),z(n)))/(1-h*fprime(z(n)));
end
y(k+1)=z(n+1);
end
plot (x,y,'r')
grid on

Si on a du y' et pas de y'' , ça devient bêtement un problème de cauchy, non ? Et hop Euler explicite, implicite ou heun ? Ca me semble logique...

Et pour le tir au but, je suis d'accord avec toi... vive les matrices qui marchent tout le temps et qui sont plus simples!

Nemo a écrit:
end
xD'

bête petite question, ça permet quoi exatement de terminé le programme par xD' (comparé à x' qui affiches toutes les valeurs de x en colonne et non en ligne) ?
merchiiii

si tu parle du xD pour la dichotomie, c'est juste que c'est le nom donné au x trouvé par dichotomi, comme yavai déja x1, x2, ca change un peu! donc x' ou xD' c'est la meme chose!

lol j'ai rien dit, j'ai mal lu le programme sryyyyyyyyyy

edit: yep j'ai vu Nemo, merci quand même

ptite question pour cette méthode des matrice (qui doit bien avoir un nom un peu plus précis que ca j'imagine...) donc jvoudrais savoir comment on discrédirte un Y', au tp il donnait pour y''=y(n+1)-2y(n)+y(n-1)/h^2 ok c'est nikel pour cet exo mais si jamais on a du y'!!? qu'est ce qui ce pass?
si quelqu'un pouvait m'renseigner ce serai bien aimable de ca part!
merci[quote]

Pour répondre à Nemo sur la question du y', tu le remplace par y(n+1) -y(n)/h. Soit la définition dune dérivée.

j'ai 2 pitites questions pour l'exo n°2
dans la méthode de matrice tridiagonale, le N est bien choisi de façon arbitraire non?
et est-ce que quelqu'un aurait la résolution de cet exo mais avec la méthode des tirs au but? parce que j'y comprend queuuuuud

merci bieng!!

yep fo choisir un N au pif, mais on dirait que plus il est grand, plus le graph est précis et "fluide"! si tu prend genre N=400 ca donne un joli graphique!!

et pour le tir au but, voila mon programme pour la seance 11, exo3. il m'a l'air juste, vu que quand je le compare a la méthode des matrices ca coîcide plus ou moins! enfin voila:

%tp11 tir au but-->exo3
clc
clear
f=inline('-sin(x)*u1+x-1')
C1=100; C2=-100;
u1(1)=1; u2(1)=C1; x(1)=0;
N=50; h=pi/N;
for n=1:N
x(n+1)=x(1)+h*n;
u1(n+1)=u1(n)+h*u2(n);
u2(n+1)=u2(n)+h*f(u1(n),x(n));
end
plot(x,u1); hold on; grid on;
T1=u1(n)+1;
pause

u2(1)=C2;
for n=1:N
x(n+1)=x(1)+h*n;
u1(n+1)=u1(n)+h*u2(n);
u2(n+1)=u2(n)+h*f(u1(n),x(n));
end
plot(x,u1,'r'); hold on; grid on;
T2=u1(n)+1;
pause

Cm=C1-((C1-C2)/(T1-T2))*T1;
u2(1)=Cm;
for n=1:N
x(n+1)=x(1)+h*n;
u1(n+1)=u1(n)+h*u2(n);
u2(n+1)=u2(n)+h*f(u1(n),x(n));
end
plot(x,u1,'y'); hold on; grid on;

jsé pas trop pq l'ordi fait un peu l'malin, il écrit en petit, donc deso pour la petitesse de la réponse!

de rien pour la petitesse mais merci pour la réponse

heu pour vos histoire de matrice, cette méthode ne marche pas a chaque fois. Cette méthode ne marche que si on a un problème linéaire, tandis que la méthode du tir au but marche pour les 2 types de problèmes (linéaire et non linéaire)