Harj. 7, LV, teht. 1 ===================== > Euler:=proc(f,a,b,ya,m) local n,h,t,y; h:=evalf((b-a)/m); t[0]:=a;y[0]:=ya; for n from 0 to m do y[n+1]:=y[n]+h*f(t[n],y[n]); t[n+1]:=t[n]+h; od; seq([t[n],y[n]],n=0..m); end: # # Proseduuri Euler haettiin www:stä Maple-koodeja-kohdasta > -------------------------------------------------------------------------------- # # Teht. 1 # ------ # Tarkoitus oli tehdä AV, teht. 2 > f:=(t,y)->-5*t^4*y^2; # Oikea Maple-funktio (ei siis pelkkä lauseke) 4 2 f := (t,y) -> - 5 t y -------------------------------------------------------------------------------- > with(DEtools): -------------------------------------------------------------------------------- > ?DEplot -------------------------------------------------------------------------------- > sk:=DEplot(diff(y(t),t)=f(t,y(t)),[t,y],t=0..2,y=-3..3,arrows=THIN,stepsize=0.1): ## Versio V.3 Pelkkä suuntakenttä ( Versiossa V.4 on kai näin:) sk:=DEplot(diff(y(t),t)=f(t,y(t)),[y],t=0..2,y=-3..3,stepsize=0.1): Tarkka analyyttinen ratkaisu: -------------------------------------------------------------------------------- > dsolve({diff(y(t),t)=-5*t^4*y(t)^2,y(0)=1},y(t)); 1 y(t) = - -------- 5 - t - 1 ----------------------------------------------------------------------------- Annetaan alkuarvoksi parametri c: > dsolve({diff(y(t),t)=-5*t^4*y(t)^2,y(0)=c},y(t)); 1 y(t) = - ---------- 5 - t - 1/c -------------------------------------------------------------------------------- > analratk:=simplify("); c analratk := y(t) = -------- 5 t c + 1 -------------------------------------------------------------------------------- > analparvi:=seq(rhs(analratk),c=-3..3): -------------------------------------------------------------------------------- > with(plots): -------------------------------------------------------------------------------- > parvi:=plot({analparvi},t=0..2,y=-4..4): -------------------------------------------------------------------------------- > display({sk,parvi}); -------------------------------------------------------------------------------- > display(parvi); -------------------------------------------------------------------------------- > m:=10:a:=0:b:=5:ya:=1: -------------------------------------------------------------------------------- > Euler(f,a,b,ya,m); [0, 1], [.5000000000, 1], [1.000000000, .8437500000], [1.500000000, -.9360351570], [2.000000000, -12.02495814], 10 [2.500000000, -5796.009693], [3.000000000, -.3280643331*10 ], 22 46 [3.500000000, -.2179430686*10 ], [4.000000000, -.1781961468*10 ], 94 190 [4.500000000, -.2032247471*10 ], [5.000000000, -.4233925845*10 ] -------------------------------------------------------------------------------- > m:=20:eu:=Euler(f,a,b,ya,m):eu[m+1]; [5.000000000, .0002694569206] -------------------------------------------------------------------------------- > m:=15:eu:=Euler(f,a,b,ya,m):eu[m+1]; 527 [4.999999996, -.8015719063*10 ] -------------------------------------------------------------------------------- > m:=16:eu:=Euler(f,a,b,ya,m):eu[m+1]; [5.000000000, .0002591634589] -------------------------------------------------------------------------------- > for m from 10 to 20 do m,Euler(f,a,b,ya,m)[m+1],evalf(subs({c=1,t=5},analratk)); od; 190 10, [5.000000000, -.4233925845*10 ], y(5) = .0003198976328 238 11, [5.000000004, -.1666179016*10 ], y(5) = .0003198976328 12, [5.000000003, .0002435761474], y(5) = .0003198976328 181 13, [5.000000004, -.2851871673*10 ], y(5) = .0003198976328 515 14, [4.999999998, -.5001295789*10 ], y(5) = .0003198976328 527 15, [4.999999996, -.8015719063*10 ], y(5) = .0003198976328 16, [5.000000000, .0002591634589], y(5) = .0003198976328 17, [4.999999999, .0002621121860], y(5) = .0003198976328 18, [5.000000003, .0002647873845], y(5) = .0003198976328 19, [5.000000004, .0002672255682], y(5) = .0003198976328 20, [5.000000000, .0002694569206], y(5) = .0003198976328 -------------------------------------------------------------------------------- > m:=100:eu:=Euler(f,a,b,ya,m):eu[m+1]; [5.000000000, .0003084149364] -------------------------------------------------------------------------------- > m:=200:eu:=Euler(f,a,b,ya,m):eu[m+1]; [5.000000000, .0003140578473] -------------------------------------------------------------------------------- > m:=1000:eu:=Euler(f,a,b,ya,m):eu[m+1]; [5.000000000, .0003187135851] -------------------------------------------------------------------------------- > Mistä johtuu 'räjähdysalttius'? ================================ No jos ajatellaan suuntakenttää, niin ihan diff. yhtälön lausekkeesta näkyy, että yksikköneliön sisäpuolella suuntakenttä kulkee miltei vaakasuorassa ja aina, jos ollaan riittävän lähellä t-akselia (tai y-akselia). Sensijaan, jos t>1 ja y>1, niin derivaatta on itseisarvpoltaan hyvin suuri (ainakin jos toinen on vähänkään iso), suuntakenttä siis miltei pystysuorassa. Jos nyt tällaisella pystysuoralla alueella lähdetään tangentin suuntaan liian pitkä askel, niin joudutaan hyvin kauas alas ja suuntakentän taipuminen vaakasuoraksi t-akselin lähellä jää kokonaan huomaamatta. Seuraavalla askeleella tilanne vain pahenee.