Universite de Nice L3MASS annee Departement de Mathematiques NOM Date PRENOM Groupe

Publié par

Universite de Nice L3MASS, annee 2011-2012 Departement de Mathematiques NOM : Date : . PRENOM : Groupe : . Feuille-question du TP 3 Les algorithmes d'Euler et de Runge-Kutta du 2eme ordre 1 Programmation de l'algorithme d'Euler On appelle algorithme de resolution d'une equation differentielle ordinaire y? = f(t, y) une fonction (t, y) 7? ?(t, y ;h) qui doit etre une bonne approximation y˜(t + h) de la solution exacte y de l'equation qui verifie y(t + h) = y. Le nombre h s'appelle le pas d'integration. L'idee est de partir d'une condition initiale (t0, y0) et de considerer la suite des ti = ti?1 + h = t0 + i ? h et des yi = ?(ti?1, yi?i, h), en esperant que yi soit une bonne approximation de la valeur de la solution y(ti) telle que y(t0) = y0. On est satisfait si, tout chaque T fixe, et pour h := T/N la limite, lorsque N devient grand (et donc le pas h devient petit), la difference ∆ = y(T ) ? yN entre valeur exacte et approximation tends vers zero, et d'autant plus si cette convergence est “rapide”, c'est a dire qu'il n'est pas necessaire de choisir N trop grand pour atteindre une precision souhaitee.

  • universite de nice l3mass

  • solution exacte

  • methode du point milieu pour l'integrale

  • methode de runge-kutta

  • feuille-question du tp

  • methode d'euler


Publié le : mardi 19 juin 2012
Lecture(s) : 28
Source : math.unice.fr
Nombre de pages : 4
Voir plus Voir moins
Universit´edeNice De´partementdeMathe´matiques NOM : PRENOM :
Date : Groupe :
L3MASS,anne´e2011-2012 . .
Feuille-question du TP 3 LesalgorithmesdEuleretdeRunge-Kuttadu2`emeordre
1 Programmationde l’algorithme d’Euler 0 Onappellealgorithmedere´solutiondunee´quationdi´erentielleordinairey=f(t, y) une fonction (t, y)7→Φ(t, y;hxiroppeaontimaodtieˆrtueenobnn)quiy˜(t+h) de la solution exacteyntaoi´equdelquiv´eriey(t+h) =y. Le nombrehs’appelle lesdpataoi.nitne´rgdepartirid´eeestLenudcondition initiale(t0, y0td)eesitedalusrereis´dcenoti=ti1+h=t0+ihet desyi= Φ(ti1, yii, h), en esp´erantqueyisoit une bonne approximation de la valeur de la solutiony(ti) telle quey(t0) =y0. On est satisfait si, tout chaqueTxrtpou´e,eh:=T /Nla limite, lorsqueNdevient grand (et donc le pas hdveeitnΔecn=adi´erepetit),ly(T)yNrtnectxataeealevreeue´zsrevste,ormaxiroppndteonti dautantplussicetteconvergenceestrapide,cest`adirequilnestpasn´ecessairedechoisirNtrop grandpouratteindreunepre´cisionsouhait´ee.Lalgorithmedeladroitetantentedˆua`Euler(1707-1783) consiste`aposer Φ(t, y, h) =y+hf(t, y).(1) 0 Nousallonslexp´erimentertoutdabordsurle´quationdie´rentielley=ydont voici le codescilab : clear ; xset("window",0) ; function ff=f(y); ff=-y; endfunction; function ff=ff(t,y); ff=f(y); endfunction; Les fonctionsf,ff, etfffsnoisvettronsdersioqe´ertonidnoitautiener´inquleelequ´eedndpe deyet pas det, mais la commandeodeexige une fonction de deux variables (t, yua;q)a`tnfff, elle t retourne le vecteur colonne(1,y) =(1,-y)’lorsquev=(t,y) ;elpeleetrmredemprpe´estnreelcah de directions” function ff=fff(t,v); ff=[1,-v(2)]’; endfunction; xr=0 :0.2 :2;yr=0 :0.4 :4;fchamp(fff,0,xr,yr) ; Saisissezleslignesci-dessuspouravoirunerepr´esentationg´eome´triquedenotree´quationdie´rentielle 0 y=yup,euedalosulitnoerenudrtnese´rp´engioatiqtr´eomlsgisiel-iedencsspoussouposerdisissuede (0,-a`-eridlleteuqe4)(cesty(0) = 4). Tmax=2 ;y0=4 ; N=100 ;petitpas=Tmax/N ; M=10 ;grandpas=Tmax/M ; t=0 :petitpas :Tmax; y=ode(y0,0,t,ff) ; plot(t,y) ;
Voici comment coder l’algorithme d’Euler sous forme d’une fonction : //Euler(t0,y0,pas) est une approximation de y(t0+pas) pour y(t0)=y0. function y=Euler(t0,y0,pas); y=y0+pas*ff(t0,y0) ; endfunction ; Saisissez-ladansscilab.Linstructionsuivantepermetalorsderepre´senterlepremierpasdelalgo-rithme, pourh= 0.urpermetzgrandpoemtnsaesoltniaer;eatcuxe´eerltsuilaulresdertsiveet-z2ovisiohc laetdonnerunschemadutrac´equevousobtenez. plot([0,grandpas],[y0, Euler(0,y0,grandpas)],’r-o’);
.
1
Lesinstructionssuivantesr´ep`etentcetalgorithmejusqua`atteindreT= 2 =Tmax. Executez-les et observezl´ecartentrelasolutionetsonapproximation. pas=grandpas ; tt=0 :pas :Tmax; yy=zeros(1+Tmax/pas) ;yy(1)=y0 ; for i=1 :1 :Tmax/pas; yy(i+1)=Euler(tt(i),yy(i),pas) ; end ; plot(tt,yy,’g->’) ; disp(y(N+1)-yy(1+Tmax/pas),’difference=’,pas,’pas=’) ;//notez l’ordre inverse Quele´cartΔtrouvez-vousentrelavaleury(2) de la solution et son approximationyM=yy(1+M) ? Δ0.2= Expliquez pourquoi on a bieny˜(T) =yy(M+1)=yy(1+Tmax/pas).
Repre´sentezavecsoinci-dessousledessinquevousobtenez;indiquezsurvotredessinquellevaleur approximative dey˜(0.4) vous lisez sur le dessin Quelle est la valeur exacte dey˜(0.epouithmlgorrlaapee´vuort)4crteh= 0.2 y˜(0.4) =
2
Reprenezcette´etudeavech= 0.02 : // et maintenant avec un petit pas pas=petitpas ; Utilisezlacommandeci-dessouspourobtenirquelenouveautrace´soitenrouge plot(tt,yy,’r-.’) ; Qu’observez-vous ? Notons Δ0.2et Δ0.02ne,sle´setracttuoianpptceestloutionexaentresol,2=rouep´echroh= 0.2 et h= 0.vzennoD.20:tatsesulosr´ Δ0.2= Δ0.02= Δ0.02/Δ0.2= Commente´voluecetteerreuraveclatailledupaschoisi?
02 Reprenezcette´etudepourl´equationdi´erentielley=ysur une nouvelle figure : ///////////////////////// on recommence avec -y^2 xset("window",1) ; Notons`anouveauΔ0.2et Δ0.02´scerastn,eelturrppanoitopee´hcorentluso,e=2teetulosnoitcaxe h= 0.2 eth= 0.ats:2.0zeovoDnnustlrse´ Δ0.2= Δ0.02= Δ0.02/Δ0.2= Commente´voluecetteerreuraveclatailledupaschoisi?
3
2Me´thodedeRunge-Kuttadu2´emeordre Cetteme´thodeestmeilleurequelame´thodedEuler:cestcequenousallonsexp´erimenterici.Elle estinspire´edelam´ethodedupointmilieupourlinte´graleetconsiste`aposer
Φ(t, y, h) =y+hf(t+h/2, y+hf(t, y)/2).
(2)
Ci-dessous son expression sous forme d’une fonctionscilab : /////////// On reprend y’=-y avec la methode de RK2 //RK(t0,y0,pas) est une approximation de y(t0+pas) pour y(t0)=y0. function y=RK(t0,y0,pas); // variables locales : k1,k2; k1=pas*ff(t0,y0) ; k2=pas*ff(t0+pas/2,y0+k1/2) ; y=y0+k2 ; endfunction ; Mise en oeuvre : xset("window",3) ; xr=0 :0.2 :2;yr=0 :0.4 :4;fchamp(fff,0,xr,yr) ;//pourune meilleur comprehension t=0 :petitpas :Tmax; y=ode(y0,0,t,ff) ; plot(t,y) ;//lasolution exacte plot([0,grandpas],[y0,RK(0,y0,grandpas)],’r-o’) ;//methodesur un seul pas pas=grandpas ;//aremplacer par petitpas pour h=0.02 tt=0 :pas :Tmax; yy=zeros(1+Tmax/pas) ;yy(1)=y0 ; for i=1 :1 :Tmax/pas; yy(i+1)=RK(tt(i),yy(i),pas) ; end ; plot(tt,yy,’g->’) ;//remplacer ’g->’ par ’r-.’ pour h=0.02 disp(y(N+1)-yy(1+Tmax/pas),’difference=’,pas,’pas=’) ;//notez l’ordre inverse Notonsa`nouveauΔ0.2et Δ0.02netr,sca´eeslt=rtse,2neioneoluteetsxactropp´echutolnaiouoper h= 0.2 eth= 0.zvosr´es02.Donneluatst: Δ0.2= Δ0.02= Δ0.02/Δ0.2= Commente´voluecetteerreuraveclatailledupaschoisi?
Notons pour finir que la commandeodech´eutilmhdecelaucalpporiussauleitorlgnaeneibesileudnetn t delasolution.Commeonsaitquelasolutionexacteestdonn´eepary(t) = 4ecalculez l’erreur par disp(4*exp(-2)-y(N+1),’difference=4exp(-2)-ode(2)’) ; Δ = Utilisez le “help” descilableedurqurmerentiglalotpseeuhmitore.s´liti
4
Soyez le premier à déposer un commentaire !

17/1000 caractères maximum.