Teht. 4
========

------------------------- tiedosto ff.m -- polun varrella ------------
function ypr=ff(t,y)
% NumaB tentti toukok. -97
ypr=(2./t).*y+t.^2.*exp(t);
-----------------------------------------------------------------------

>> clear t y                                                
>> t0=1;y0=0;h=0.25;y(1)=y0;t(1)=t0;n=1;                    
>>  y(n+1)=y(n)+h*ff(t(n),y(n));t(n+1)=t(n)+h;n=n+1;[t',y'] 

ans =

    1.0000         0
    1.2500    0.6796

>>  y(n+1)=y(n)+h*ff(t(n),y(n));t(n+1)=t(n)+h;n=n+1;[t',y']

ans =

    1.0000         0
    1.2500    0.6796
    1.5000    2.3148

>>  y(n+1)=y(n)+h*ff(t(n),y(n));t(n+1)=t(n)+h;n=n+1;[t',y']

ans =

    1.0000         0
    1.2500    0.6796
    1.5000    2.3148
    1.7500    5.6074

>>  y(n+1)=y(n)+h*ff(t(n),y(n));t(n+1)=t(n)+h;n=n+1;[t',y']

ans =

    1.0000         0
    1.2500    0.6796
    1.5000    2.3148
    1.7500    5.6074
    2.0000   11.6153

>> plot(t,y)
%
% Tarkka ratkaisu:
%  Kuva:
>> tt=linspace(1,2);yy=(tt.^2).*(exp(tt)-exp(1));
>> plot(tt,yy,'g')

% Lasketaan tarkka edellisissä laskentapisteissä:
>> ytarkka=(t.^2).*(exp(t)-exp(1)); 
% Virheet:
>> absv=ytarkka-y
absv =

         0    0.5268    1.6529    3.6914    7.0678    <- VASTAUS
              =====================================

>> suhtv=absv./ytarkka
Warning: Divide by zero

suhtv =

       NaN    0.4367    0.4166    0.3970    0.3783     <- VASTAUS
              =====================================

% Alkuarvopisteessä arvo = 0, joten suht. virhe = NaN, no siihen
% ei mielenkiintomme kohdistu, sehän oli annettu.

-----------------------------------------------------------------

% Jos halutaan katsoa, mitä tapahtuu pienemmällä askeleella, kannattaa
% nn:n sijasta käyttää for- tai while- rakennetta:
%
>> t025=t;y025=y;  % Talletetaan askeleella h=0.25 lasketut.
>> clear t y                                                
>> t0=1;y0=0;h=0.1;y(1)=y0;t(1)=t0;n=1; % Lasketaan askeleella h=0.1
%
>>  while t<2 , y(n+1)=y(n)+h*ff(t(n),y(n));t(n+1)=t(n)+h;n=n+1;[t',y']; 
end
>> plot(t025,y025,t,y,'g')