%-- 29.11.2012 8:21 --% help ode45 [t,y]=ode45(@vdp1,[0 20],[2 0]); plot(t,y(:,1)); odeset h=odeset help ode45 [t,y]=ode45('t1',[0,20],1); plot(t,y) hold on [t,y]=ode113('yp',[t0,tf],y0); [t,y]=ode113('t1',[0,20],1); plot(t,y,'r') [t,y]=ode15s('t1',[0,20],1); plot(t,y,'k') [t,y]=ode15s('t1b',[0,20],1); plot(t,y,'k') [t,y]=ode45t1b',[0,20],1); [t,y]=ode45('t1b',[0,20],1); plot(t,y,'k') odeset h=odeset h=set('AbsTol',1) set(h,'AbsTol',1) %Kokeillaan vähän odeset(h,'AbsTol',1) odeset(h,'AbsTol',true) odeset('RelTol', 1e-20, 'AbsTol', 1e-4) odeset('RelTol', 1e-20, 'AbsTol', 1e-20) [t,y]=ode45('t1b',[0,20],1); plot(t,y,'k') [t,y]=ode45('t1',[0,20],1); h=odeset('RelTol', 1e-20, 'AbsTol', 1e-20) [t,y]=ode45('t1',[0,20],1,h); h=odeset('RelTol', 1e-10, 'AbsTol', 1e-10) [t,y]=ode45('t1',[0,20],1,h); plot(t,y,'k') h=odeset('RelTol', 1, 'AbsTol', 1) [t,y]=ode45('t1',[0,20],1,h); h=odeset('RelTol', 1e-1, 'AbsTol', 1e-1) [t,y]=ode45('t1',[0,20],1,h); plot(t,y,'k') [t,y]=ode113('t1',[0,20],1,h); hold on plot(t,y,'r') [t,y]=15s('t1',[0,20],1,h); [t,y]=ode15s('t1',[0,20],1,h); plot(t,y,'b') h=odeset('RelTol', 1e-.5, 'AbsTol', 1e-.5) h=odeset('RelTol', 1e-0.5, 'AbsTol', 1e-0.5) h=odeset('RelTol', 1e-1, 'AbsTol', 1e-15) h=odeset('RelTol', 1e-1, 'AbsTol', 1e-1) h=odeset('RelTol', 1e-1, 'AbsTol', 1e-15) [t,y]=ode15s('t1',[0,20],1,h); plot(t,y,'b') h=odeset('RelTol', 1e-15, 'AbsTol', 1e-1) [t,y]=ode15s('t1',[0,20],1,h); h=odeset('RelTol', 1e-10, 'AbsTol', 1e-1) [t,y]=ode15s('t1',[0,20],1,h); plot(t,y,'b') clear %Numeerinen analyysi s2012, Laskuharjoitus 12 888888888888 88 88 88 88 ,d ,d88 88 88 88 888888 88 ,adPPYba, 88,dPPYba, MM88MMM 88 88 a8P_____88 88P' "8a 88 88 88 8PP""""""" 88 88 88 88 88 "8b, ,aa 88 88 88, 888 88 88 `"Ybbd8"' 88 88 "Y888 888 88 %Difyhtälön funktio on tallennettu kansioon nimellä t1.m. Funktio ilman sinitermiä on tallennettu nimellä t1b.m. %Kokeillaan aluksi ratkaista yhtälö eri menetelmillä. [t,y]=ode45('t1',[0,20],1); plot(t,y) a=t;y=b; a=t;b=y; [t,y]=ode113('t1',[0,20],1); hold on plot(t,y-1,'k') c=t;d=y; [t,y]=ode15s('t1',[0,20],1); plot(t,y-2,'b') plot(t,y-2,'m') title('Harjoitus 12, tehtävä 1, ratkaisuja eri menetelmillä') legend('ode45','ode113 -1','ode15s -2') xtitle('t') xlabel('t') ylabel('x') print -djpeg 't1.jpg' %Piirretään samat ratkaisut ilman sinitermiä [t,y]=ode45('t1b',[0,20],1); plot(t,y) hold on [t,y]=ode113('t1b',[0,20],1); plot(t,y-1,'k') [t,y]=ode15s('t1b',[0,20],1); plot(t,y-2,'m') title('Harjoitus 12, tehtävä 1, ratkaisuja eri menetelmillä, kun sinitermiä ei ole') legend('ode45','ode113 -1','ode15s -2') xlabel('t') ylabel('x') print -djpeg 't1b.jpg' %Kokeillaan vielä erilaisia toleransseja odeset %tässä on kai lista ominaisuuksista h=odeset %tässä on erilainen lista AbsTol odeset('RelTol', 1e-20, 'AbsTol', 1e-20) h=odeset('RelTol', 1e-20, 'AbsTol', 1e-20) [t,y]=ode45('t1',[0,20],h); clear [t,y]=ode45('t1',[0,20],h); h=odeset('RelTol', 1e-20, 'AbsTol', 1e-20) [t,y]=ode45('t1',[0,20],h); [t,y]=ode45('t1',[0,20],1,h); %Matlab ilmoittaa nostaneensa RelTolia, joka oli Matlabille liian pieni plot(t,y) %Tulee hyvä ratkaisu h=odeset('RelTol', 1, 'AbsTol', 1); [t,y]=ode45('t1',[0,20],1,h); %Nyt kun laitettiin liian iso toleranssi, Matlab kertoo, että ratkaiseminen ei onnistunut. %Kokeillaan vielä, kumpi toleranssi on merkityksellisempi h=odeset('RelTol', 1e-10, 'AbsTol', 1); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) %ratkaisualgoritmi menee läpi, mutta "ratkaisu" on jotain ihan kauheata print -djpeg 't1_kauhea.jpg' h=odeset('RelTol', 1, 'AbsTol', 1e-10); [t,y]=ode45('t1',[0,20],1,h); %RelTol oli iso ja AbsTol oli normaali, mutta laskenta ei mennyt läpi %Siis Matlab ei voi käsitellä kovin pieniä toleransseja %Jos RelTol on iso, laskenta menee jossakin vaiheessa pieleen %Jos RelTol on normaali ja AbsTol iso, tulee epätarkkoja ratkaisuja %Aiheesta voisi tehdä videon %Kasvatetaan videolla AbsTolia normaalista isoksi for j=1:36 h=odeset('RelTol', 10^(-10+j/4), 'AbsTol', 10^(-10+j/4)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end movie(F) %tuli ihan hauska video, mutta liian lyhyt save 'video.mat' [t,y]=ode113('t1b',[0,20],1); plot(t,y-1,'k') [t,y]=ode15s('t1b',[0,20],1); plot(t,y-2,'m') title('Harjoitus 12, tehtävä 1, ratkaisuja eri menetelmillä, kun sinitermiä ei ole') legend('ode45','ode113 -1','ode15s -2') xlabel('t') ylabel('x') print -djpeg 't1b.jpg' for j=1:370 h=odeset('RelTol', 10^(-10+j/40), 'AbsTol', 10^(-10+j/4)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end for j=1:370 h=odeset('RelTol', 10^(-10+j/40), 'AbsTol', 10^(-10+j/40)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end clear for j=1:370 h=odeset('RelTol', 10^(-10+j/40), 'AbsTol', 10^(-10+j/40)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end save 'video2.mat' clear load video.mat movie(F,1,1) clear for j=1700:1800 h=odeset('RelTol', 10^(-10+j/200), 'AbsTol', 10^(-10+j/200)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end movie(F) clear for j=1600:1750 h=odeset('RelTol', 10^(-10+j/200), 'AbsTol', 10^(-10+j/200)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end clear for j=1500:1700 h=odeset('RelTol', 10^(-10+j/200), 'AbsTol', 10^(-10+j/200)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end save 'video3.mat' for j=1700:1800 h=odeset('RelTol', 10^(-10), 'AbsTol', 10^(-10+j/200)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end for j=1700:1900 h=odeset('RelTol', 10^(-10), 'AbsTol', 10^(-10+j/200)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end movie(F) clear for j=1700:1850 h=odeset('RelTol', 10^(-10), 'AbsTol', 10^(-10+j/200)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end clear for j=1700:1800 h=odeset('RelTol', 10^(-10), 'AbsTol', 10^(-10+j/200)); [t,y]=ode45('t1',[0,20],1,h); plot(t,y) hold on [t,y]=ode113('t1',[0,20],1,h); plot(t,y-1,'r') [t,y]=ode15s('t1',[0,20],1,h); plot(t,y-2,'k') F(j)=getframe; close end save 'video4.mat' clear 888888888888 88 ad888888b, 88 88 ,d d8" "88 88 88 88 a8P 88 ,adPPYba, 88,dPPYba, MM88MMM ,d8P" 88 a8P_____88 88P' "8a 88 a8P" 88 8PP""""""" 88 88 88 a8P' 88 "8b, ,aa 88 88 88, 888 d8" 88 `"Ybbd8"' 88 88 "Y888 888 88888888888 [t,y]=ode15s('t2',[0,0.99],1); [t,y]=ode45('t2',[0,0.99],1); plot(t,y) hold on [t,y]=ode113('t2',[0,0.99],1); plot(t,y+1) plot(t,y+1,'w') [t,y]=ode45('t2',[0,0.99],1); close [t,y]=ode45('t2',[0,0.99],1); plot(t,y+1) hold on [t,y]=ode113('t2',[0,0.99],1); close [t,y]=ode45('t2',[0,0.99],1); plot(t,y) hold on [t,y]=ode113('t2',[0,0.99],1); plot(t,y+10) plot(t,y+10,'r') [t,y]=ode15s('t2',[0,0.99],1); plot(t,y+20,'k') title('Harjoitus 12, tehtävä 3, ratkaisuja eri menetelmillä') legend('ode45','ode113 +10','ode15s +20') xlabel('t') ylabel('x') print -djpeg 't2.jpg' title('Harjoitus 12, tehtävä 2, ratkaisuja eri menetelmillä') print -djpeg 't2.jpg' close [t,y]=ode45('t2',[0,2],1); [t,y]=ode113('t2',[0,2],1); [t,y]=ode15s('t2',[0,2],1); %mikään menetelmä ei mene arvoista t<1 arvoihin t>1 %Menetelmät näyttävät käyttäytyvän samalla tavalla. 888888888888 88 ad888888b, 88 88 ,d d8" "88 88 88 88 a8P 88 ,adPPYba, 88,dPPYba, MM88MMM aad8" 88 a8P_____88 88P' "8a 88 ""Y8, 88 8PP""""""" 88 88 88 "8b 88 "8b, ,aa 88 88 88, 888 Y8, a88 88 `"Ybbd8"' 88 88 "Y888 888 "Y888888P' h=odeset('outputfcn',@odephas2); clear t clear y [t,y]=ode15s('t3',[0,2],[1,1]); plot(t,y) [t,y]=ode15s('t3',[0,2],[1,1],h); hold on [t,y]=ode15s('t3',[0,2],[2,2],h); [t,y]=ode15s('t3',[0,2],[0,2],h); [t,y]=ode15s('t3',[0,2],[-1,-1],h); [t,y]=ode45('t3',[0,2],[0.1,0.1],h); [t,y]=ode45('t3',[0,10],[0.1,0.1],h); close [t,y]=ode45('t3',[0,10],[0.1,0.1],h); [t,y]=ode45('t3',[0,10],[0.1,1.1],h); [t,y]=ode45('t3',[0,1],[0.1,1.1],h); [t,y]=ode45('t3',[0,100],[0.1,1.1],h); [t,y]=ode45('t3',[0,100],[2.1,2.1],h); [t,y]=ode45('t3',[0,100],[2.0,2.1],h); [t,y]=ode45('t3',[0,100],[2.1,2.1],h); [t,y]=ode45('t3',[0,100],[2.1,1.9],h); [t,y]=ode45('t3',[0,100],[2.1,1.8],h); [t,y]=ode45('t3',[0,10],[0.1,1.1],h); [t,y]=ode45('t3',[0,10],[0.1,0.1],h); hold on [t,y]=ode45('t3',[0,10],[-1.1,0.99],h); [t,y]=ode45('t3',[0,10],[-1.1,1.99],h); [t,y]=ode45('t3',[0,10],[-1.2,3],h); clear [t,y]=ode45('t3',[0,10],[-1.2,3],h); h=odeset('outputfcn',@odephas2); [t,y]=ode45('t3',[0,10],[-1.2,3],h); close [t,y]=ode45('t3',[0,10],[-1.2,3],h); [t,y]=ode45('t3',[0,10],[-1.1,1.99],h); [t,y]=ode45('t3',[0,10],[-1.1,3],h); [t,y]=ode45('t3',[0,10],[-1.1,2.5],h); [t,y]=ode45('t3',[0,10],[-1.1,2.4],h); [t,y]=ode45('t3',[0,10],[-1.1,2.3],h); [t,y]=ode45('t3',[0,10],[-1.1,2.2],h); [t,y]=ode45('t3',[0,10],[-1.1,2.1],h); [t,y]=ode45('t3',[0,10],[-1.1,2.05],h); [t,y]=ode45('t3',[0,10],[-1.1,2.09],h); [t,y]=ode45('t3',[0,10],[-1.1,2.08],h); [t,y]=ode45('t3',[0,10],[-1.1,2.07],h); [t,y]=ode45('t3',[0,10],[-1.1,2.06],h); [t,y]=ode45('t3',[0,10],[-1.1,2.065],h); [t,y]=ode45('t3',[0,10],[-1.1,2.064],h); [t,y]=ode45('t3',[0,10],[-1.1,2.063],h); [t,y]=ode45('t3',[0,10],[-1.2,2.063],h); [t,y]=ode45('t3',[0,10],[-1.2,1.8],h); [t,y]=ode45('t3',[0,10],[-1.1,1.8],h); [t,y]=ode45('t3',[0,10],[-1.1,2.5],h); [t,y]=ode45('t3',[0,10],[-0.9,2.5],h); [t,y]=ode45('t3',[0,10],[-0.9,1.5],h); [t,y]=ode45('t3',[0,10],[-0.9,1.4],h); [t,y]=ode45('t3',[0,10],[-0.9,1.8],h); [t,y]=ode45('t3',[0,10],[-1.1,2.5],h); [t,y]=ode45('t3',[0,10],[-0.1,1.1],h); [t,y]=ode45('t3',[0,10],[0.1,1.1],h); [t,y]=ode45('t3',[0,10],[0.1,0.9],h); [t,y]=ode45('t3',[0,10],[0.1,0.8],h); [t,y]=ode45('t3',[0,10],[-10,10],h); [t,y]=ode45('t3',[0,10],[10,10],h); [t,y]=ode45('t3',[0,10],[10,-10],h); [t,y]=ode45('t3',[0,10],[-10,-10],h); [t,y]=ode45('t3',[0,10],[10,-10],h); [t,y]=ode45('t3',[0,10],[-10,-10],h); [t,y]=ode45('t3',[0,10],[-10,10],h); [t,y]=ode45('t3',[0,10],[10,10],h); [t,y]=ode45('t3',[0,10],[-10,-10],h); hold on [t,y]=ode45('t3',[0,10],[10,10],h); plot(t,y) [t,y]=ode45('t3',[0,10],[-10,-10],h); hold on [t,y]=ode45('t3',[0,10],[10,10],h); [t,y]=ode45('t3',[0,10],[-10,10],h); [t,y]=ode45('t3',[0,10],[10,10],h); [t,y]=ode45('t3',[0,10],[-10,10],h); [t,y]=ode45('t3',[0,10],[-10,-10],h); [t,y]=ode45('t3',[0,10],[10,-10],h); axis=([-15 15 -15 15]) axis([-15 15 -15 15]) axes([-15 15 -15 15]) axis([-15 15 -15 15]) print -djpeg 't3.jpg' title('Harjoitus 12, tehtävä 3, alkuarvot [\pm 10,\pm 10]') title('Harjoitus 12, tehtävä 3, alkuarvot [\pm 10,\pm 10], t=0..10') print -djpeg 't3.jpg' close [t,y]=ode45('t3',[0,10],[-1.1,2.5],h); hold on [t,y]=ode45('t3',[0,10],[-0.9,2.5],h); [t,y]=ode45('t3',[0,10],[-0.9,1.5],h); [t,y]=ode45('t3',[0,10],[-0.9,1.4],h); [t,y]=ode45('t3',[0,10],[-0.9,1.8],h); [t,y]=ode45('t3',[0,10],[-1.1,2.5],h); [t,y]=ode45('t3',[0,10],[-0.1,1.1],h); [t,y]=ode45('t3',[0,10],[0.1,1.1],h); [t,y]=ode45('t3',[0,10],[0.1,0.9],h); [t,y]=ode45('t3',[0,10],[0.1,0.8],h); [t,y]=ode45('t3',[0,10],[-1.1,2.5],h); hold on [t,y]=ode45('t3',[0,10],[-0.9,2.5],h); [t,y]=ode45('t3',[0,10],[-0.9,1.5],h); [t,y]=ode45('t3',[0,10],[-0.9,1.4],h); [t,y]=ode45('t3',[0,10],[-0.9,1.8],h); [t,y]=ode45('t3',[0,10],[-1.1,2.5],h); [t,y]=ode45('t3',[0,10],[-0.1,1.1],h); [t,y]=ode45('t3',[0,10],[0.1,1.1],h); [t,y]=ode45('t3',[0,10],[0.1,0.9],h); [t,y]=ode45('t3',[0,10],[0.1,0.8],h); axis([-15 15 -15 15]) [t,y]=ode45('t3',[0,10],[-1.1,2.5],h); hold on [t,y]=ode45('t3',[0,10],[-0.9,2.5],h); [t,y]=ode45('t3',[0,10],[-0.9,1.5],h); [t,y]=ode45('t3',[0,10],[-0.9,1.4],h); [t,y]=ode45('t3',[0,10],[-0.9,1.8],h); [t,y]=ode45('t3',[0,10],[-1.1,2.5],h); [t,y]=ode45('t3',[0,10],[-0.1,1.1],h); [t,y]=ode45('t3',[0,10],[0.1,1.1],h); [t,y]=ode45('t3',[0,10],[0.1,0.9],h); [t,y]=ode45('t3',[0,10],[0.1,0.8],h); axis([-15 15 -15 15]) axis equal clear axis axis([-15 15 -15 15]) axis([-5 5 -5 5]) print close [t,y]=ode113('t3',[0,10],[0.1,0.8],h); [t,y]=ode113('t3',[0,10],[-10,-10],h); [t,y]=ode15s('t3',[0,10],[-10,-10],h); [t,y]=ode15s('t3',[0,10],[-1,-1],h); for j=1:20 a=randn(1,2); [t,y]=ode15s('t3',[0,10],a,h); end %Valaiseva silmukka %Tuntuu, että kun alkuarvoja arpoo normaalijakautuneesti, saa heti näkyviin sellaisia kohtia, joita ei olisi muuten keksinyt. %Kaikki menetelmät toimivat ongelmitta 888888888888 88 ,d8 88 88 ,d ,d888 88 88 88 ,d8" 88 88 ,adPPYba, 88,dPPYba, MM88MMM ,d8" 88 88 a8P_____88 88P' "8a 88 ,d8" 88 88 8PP""""""" 88 88 88 8888888888888 88 "8b, ,aa 88 88 88, 888 88 88 `"Ybbd8"' 88 88 "Y888 888 88 [t,y]=ode15s('t3',[0,10],[1,1]); plot(t,y) h=odeset('outputfcn',@odephas2); [t,y]=ode45('t3',[0,10],[1,1]); [t,y]=ode45('t3',[0,10],[1,1],h); [t,y]=ode45('t3',[0,10],[1,0],h); [t,y]=ode45('t3',[0,10],[1,1],h); [t,y]=ode45('t3',[0,10],[0,0],h); [t,y]=ode45('t3',[0,10],[0.1,0.1],h); [t,y]=ode45('t4',[0,10],[0.1,0.1],h); %Arvotaan lukuja %Arvotaan alkuarvoja for j=1:20 a=randn(1,2); [t,y]=ode45('t4',[0,10],a,h); end %Ratkaisut kiertävät kehää for j=1:20 a=randn(1,2); [t,y]=ode113('t4',[0,10],a,h); end %Nyt ratkaisut näyttävät vähemmän sileiltä, mutta muuten samanlaisilta for j=1:20 a=randn(1,2); [t,y]=ode15s('t4',[0,10],a,h); end %Samanlaisia ratkaisuja näkyy taaskin x=linspace(1,10); plot(x,t) x=linspace(1,10,length(y(1))) x=linspace(1,10,length(y'(1))); size(y) [c,d]=size(y); c [c,d]=size(y); x=linspace(1,10,length(c)); plot(x,y) plot(x,y(1)) plot(x,y(:,1)) y x x=linspace(0,10,length(c)); x x=linspace(0,10,c); plot(x,y(:,1)) [t,y]=ode15s('t4',[0,100],a,h); [c,d]=size(y); x=linspace(0,10,c); plot(x,y(:,1)) x=1:c; plot(x,y(:,1)) [t,y]=ode15s('t4',[0,30],a,h); [c,d]=size(y); x=linspace(0,10,c); plot(x,y(:,1)) title('Harjoitus 12, tehtävä 4, x1 näyttää Airyn funktiolta') print -djpeg 't4.jpg' plot(x,y(:,2)) title('Harjoitus 12, tehtävä 4, x2 värähtelee') print -djpeg 't4b.jpg' airy(0) -airy(2,0) t=linspace(-10,10); y=airy(-t); plot(t,y) t=linspace(0,10); y=airy(-t); plot(t,y) t=linspace(-10,0); y=airy(-t); plot(t,y) y=real(airy(-t)); plot(t,y) t=linspace(1,2); plot(t,y) y=real(airy(-2,t)); y=airy(2,-t); plot(t,y) t=linspace(-10,-1); y=airy(2,-t); plot(t,y) t=linspace(-10,-1); y=airy(-t); plot(t,y) i^2 real(1+2*i) abs(1+2*i) t=linspace(-10,-1); y=airy(-t); x=abs(y) x=abs(y); plot(t,x) t=linspace(-10,-5); y=airy(-t); x=abs(y); plot(t,x) t=linspace(-1,0); y=airy(-t); plot(t,x) t=linspace(0,1); y=airy(-t); x=abs(y); plot(t,x) t=linspace(0,10); y=airy(-t); x=abs(y); plot(t,x) hold on y=airy(2,-t); x=abs(y); plot(t,x) plot(t,x,'r') title('Harjoitus 12, tehtävä 4, Airyn funktioita') legend('airy(-t)','airy(2,-t)') legend('airy(-t)','','airy(2,-t)') print -djpeg 't4c.jpg' [t,y]=ode45('t4',[0,30],[airy(0),-airy(2,0)],h); close [t,y]=ode45('t4',[0,30],[airy(0),-airy(2,0)],h); x1=y(:,1); x2=y(:,2); t=linspace(0,30,557) t=linspace(0,30,557); plot(t,x,'b') t=linspace(0,30,557)'; plot(t,x1,'b') plot(t,x2,'r') hold on plot(t,x1,'b') y=airy(-t); plot(t,y+1,'k') z=airy(2,-t); plot(t,z+1,'m') y2=abs(airy(-t)); z2=abs(airy(2,-t)); plot(t,y2-5,'g') plot(t,z2-5,'p') plot(t,z2-5,'wp') plot(t,y2-5,'g') plot(t,z2-5,'k') title('Harjoitus 12, tehtävä 4, Airyn funktioita') legend('airy(-t) numeerisesti ratkaistuna','','airy(2,-t) numeerisesti ratkaistuna') legend('airy(-t) numeerisesti ratkaistuna','airy(2,-t) numeerisesti ratkaistuna') print -djpeg 't4d.jpg' 888888888888 88 8888888888 88 88 ,d 88 88 88 88 88 ____ 88 ,adPPYba, 88,dPPYba, MM88MMM 88a8PPPP8b, 88 a8P_____88 88P' "8a 88 PP" `8b 88 8PP""""""" 88 88 88 d8 88 "8b, ,aa 88 88 88, 888 Y8a a8P 88 `"Ybbd8"' 88 88 "Y888 888 "Y88888P" clear h=odeset('outputfcn',@odephas2); h=odeset('outputfcn',@odephas2,'nonnegative',[1 2 3]); [t,y]=ode45('t5',[0,30],[1 1 1]); %oli hidas [t,y]=ode113('t5',[0,30],[1 1 1]); %oli hidas [t,y]=ode15s('t5',[0,30],[1 1 1]); %valmistui heti y1=y(:,1); y2=y(:,2); y3=y(:,3); plot(t,y1) plot(t,y2) plot(t,y3) axis([-30 30 -30 30]) [t,y]=ode113('t5',[0,30],[1 1 1]); [t,y]=ode15s('t5',[0,30],[1 1 1],h); [t,y]=ode15s('t5',[0,30],randn(1,3),h); [t,y]=ode15s('t5',[0,30],abs(randn(1,3)),h);