%% http://www.math.hut.fi/teaching/osdylask/matlab/collection.mlb
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This is a collection of Matlab codes developed or found especially
%%%   for the needs of the "OSDYLASK 99" course
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Direction field m-file (script) 
% Source: Bäckström: Praktisk matematik med Matlab
clear variables; hold off;
xmax=4; ymax=3; delx=xmax/30; dely=ymax/30;   % *** problem dependent
[X,Y]=meshgrid(0:delx:xmax,0:dely:ymax);
  DX=ones(size(X));  % Unit x-xomponents of arrows
  DY=fode1(X,Y);      % Slopes, *** Change the name 'fode' to your m-file name
  L=sqrt(DX.^2+DY.^2); % Initial lengths
  DX1=DX./L; DY1=DY./L; % Unit length for all arrows
%figure(1); quiver(X,Y,DX1,DY1,'.'); grid on;  % The period ('.') suppresses
                                               % arrowheads
figure(1); quiver(X,Y,DX1,DY1); grid on;       % Include arrowheads
  axis([0 xmax 0 ymax]); title('Direction field');

%%% Example session (remove first %-sign)
%m=2; hold on;
%% cut/paste into Matlab-session: "up-arrow-iterate"
[T,Y]=meuler('fode1',0,3,2,m); plot(T,Y,'r',T,Y,'o');[T' Y'],m=2*m;

[t0,y0]=ginput(1);[T,Y]=meuler('fode1',t0,4,y0,m); plot(T,Y,'r');

%%%%%%% Example ode-definition files %%%%%%%
function yp=fode1(X,Y);        % file fode1.m
yp=-Y;
%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%% From this on, old y3.m-file %%%%%%%%%%%%%%%%%%%%
%==============================================
%
%
%  INTERPOLOINTI 
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lagrangen interpolaatiopolynomi
function y=lagrange(xdata,xev,j)
% Lagrangen kertojafunktio ("Cardinal function")  HA maalisk. -95
% y=lagrange(xdata,xeval,j); 
% Lasketaan xdata:aan liittyv{ j:s l-funktio pisteiss{ xeval 
% esim: y=lagrange([1 2 3 4],[1 2 3 4],0) antaa 1 0 0 0
% 
j=j+1; % Matlab-vektorissa indeksi alkaa 1:st{
y=1; tuloind=1:length(xdata);ind=tuloind~=j;tuloind=tuloind(ind);
for i=tuloind;y=y.*(xev-xdata(i))./(xdata(j)-xdata(i));end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% NUMEERINEN INTEGROINTI
%
function [S,R]=simpson(f,a,b,n,M4)
% Sy|tteet:f   - integroitava funktio
%          a,b - v{lin p{{tepisteet
%            n - "pikkuv{lien lkm. parillinen, M4: f'''':n itseisarvon. max
% Tulokset: S Simpsonarvo, R virhearvio
h=(b-a)/n;
w=2*ones(1,n-1);w(1:2:n-1)=4*ones(1,n/2); % jos n pariton, niin ruma virheilmo
w=[1,w,1];
x=a:h:b;y=feval(f,x);
S=(h/3)*w*y';
if nargout==2
  R=(M4*(b-a)*h^4)/180;
end;
%
% DIFFERENTIAALIYHTÄLÖSYSTEEMIT
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lineaarinen diffyhtälösysteemi expm-funktiota käyttäen
%-------------------------------------------------------------------
%Leikkaa ja liimaa, sijoita tiedostoon linsys.m
%--------------------------------------------------------------------
function [t,y] = linsys(A,x0,T,siz)
% Esim: A=[1 0 0;1 3 0;0 1 1];x0=[1;2;3];
%       [t,y]=linsys(A,x0,2);
%       plot(t,y)                       % aikariippuvuus
%       plot(y(:,1),y(:,2))             % faasitaso (1,2)
%       plot3(y(:,1),y(:,2),y(:,3))     % faasiavaruus (3-d)
if nargin < 4, siz = 100; end; 
t = linspace(0,T,siz);
m=size(A);m=m(1);
y = zeros(m,siz);
for i=1:siz,
y(:,i) = expm(A * t(i)) * x0;   % Tässä on asian ydin
end;
y=y';

%%%%%%%% Tavallisten diff. yhtälöiden numeriikkaa: %%%%%%%%%%%%%%
%%%%%%%%%%%%%% Eulerin menetelmä %%%%%%%%%%%%%%
%
%Nuolinäppäin-Euler (eli 'laskinmoodi')

% Esimerkki:

Yhtälö: y'=t-y^2, y(0)=1 

Kirjoita tiedostoon feu.m  asianmukainen funktiokoodi:
---------------------
function ypr=feu(t,y)
ypr=t-y.^2;
--------------------
%% Nyt voit laskea:
>> clear t y                              % Poistetaan varmuuden vuoksi.
>> t0=0;y0=0;h=0.1;y(1)=y0;t(1)=t0;n=1;   % Alkuasetukset
>> y(n+1)=y(n)+h*feu(t(n),y(n));t(n+1)=t(n)+h;n=n+1;[t',y'] % Iteraatio,toista nuolella
%% Muutaman toiston jälkeen voit tehdä:
>> plot(t,y)
%% Tästä voit jatkaa eteenpäin aika-akselilla. 
%% Jos haluat tarkentaa, aloita alkuasetusriviltä ja valitse pienempi h.
%%
% Yllä oleva nn-leikki voidaan korvata for-tai while-silmukalla. 
>>  while t<1 , y(n+1)=y(n)+h*feu(t(n),y(n));t(n+1)=t(n)+h;n=n+1;end
>> plot(t,y)
%
% Seuraava vaihe
% on kirjoittaa koko homma omaksi Matlab-funktioksi, nimeltään meuler:
%
%-------------------------------------------------------------------
%Leikkaa ja liimaa, sijoita tiedostoon meuler.m   
% (euler on varattu))
%--------------------------------------------------------------------
%
function [T,Y] = meuler(f,a,b,ya,m)
% [T,Y] = meuler(f,a,b,ya,m)
% Euler's solution for y' = f(t,y) with y(a) = ya.
% f  is the function, input.
% a  is the left endpoint, input.
% b  is the right endpoint, input.
% ya is the initial condition, input.
% m  is the number of steps, input.
% T  is the vector of abscissas, output.
% Y  is the vector of ordinates, output.
h = (b - a)/m;
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m,
  Y(j+1) = Y(j) + h*feval(f,T(j),Y(j));
  T(j+1) = a + h*j;
end

% Käyttöesimerkki:

Yhtälö: y'=t-y^2, y(0)=1 

Voit toimia seuraavan ehdotelman mukaisesti:

Kirjoita tiedostoon feu.m  asianmukainen funktiokoodi (ellet jo kirjoittanut):
---------------------
function ypr=feu(t,y)
ypr=t-y.^2;
---------------------

Matlab-istunto:

>> hold on,grid,m=3;
>> [T,Y]=meuler('feu',0,5,1,m);plot(T,Y,T,Y,'o');[T',Y'],m=2*m
>>
Toista tätä riviä (nuolinäppäimellä) ja muuta plot-komento muotoon
>> plot(T,Y)
kun m on tullut kyllin suureksi.

%=========================================================
%  Runge - Kutta 4
%=========================================================
% Leikkaa ja liimaa: sijoita tiedostoon rk4.m
%---------------------------------------------------
function [T,Y] = rk4(f,a,b,ya,m)
% [T,Y] = rk4(f,a,b,ya,m)
% Runge-Kutta solution for y' = f(t,y) with y(a) = ya.
% f  is the function, input.
% a  is the left endpoint, input.
% b  is the right endpoint, input.
% ya is the initial condition, input.
% m  is the number of steps, input.
% T  is the vector of abscissas, output.
% Y  is the vector of ordinates, output.
h = (b - a)/m;
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m,
  tj = T(j);
  yj = Y(j);
  k1 = h*feval(f,tj,yj);
  k2 = h*feval(f,tj+h/2,yj+k1/2);
  k3 = h*feval(f,tj+h/2,yj+k2/2);
  k4 = h*feval(f,tj+h,yj+k3);
  Y(j+1) = yj + (k1 + 2*k2 + 2*k3 + k4)/6;
  T(j+1) = a + h*j;
end

%----------------------------------------------------
% Käyttöesimerkit: Kuten Euler

%%%%%%% Osittaisdifferentiaaliyhtälöitä %%%%%%
%
%% Laplacen yhtälö %%%

#Tässä nyt välissä Maple-koodi

lapmatr:=proc(p)
local B,C,v;
B:=band([1,-4,1],p);
C:=diag(seq(B,i=1..p));
v:=vector(2*p+1,0);
v[1]:=1;v[2*p+1]:=1;
evalm(C+band(v,p^2))
end:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lämpöyhtälön analyyttinen ratkaisu muuttujien erottelulla.
% Written by HA s. 1995 
% Sijoita tiedostoon osdysin.m
%
function [Z,x,t]=osdysin(l,al,b,Tf,M,N,nsum)
% Lasketaan l{mp|yht{l|n sinisarjaratk. korkeusmatriisina
% Input: l=jakso,al=alpha, b=sinisarjan kerroinvektori (nollat mukana)
%        Tf=loppuaika, M=x-diskrertointi, N=t-diskr.nsum=mihin summataan
% Output: [Z,x,t] -  korkeusmatriisi, x-vektori, t-vektori
%
% Esim:
%l=10;al=sqrt(1.14);M=6;N=7;nsum=10;c=400/pi;
%clear b;nr=10;b=1./(1:nr);b=rem((1:nr),2).*b;b=c*b;Tf=10; % too fancy
% Huom: Kerroinvektori b on kätevää tehdä Maplella ja sitten vain
% leikata Matlab-istuntoon. Sarja on joka tapauksessa syytä esittää 
% muodossa,
% jossa summausindeksi juoksee kaikki luonnoll. luvut (ei siis 1,3,5,...)
% [Z,x,t]=osdysin(l,al,b,Tf,6,7,10);
%surfc(x,t,Z);
%plot(x,Z)
%plot(t,Z')
x=linspace(0,l,M);t=linspace(0,Tf,N);
Z=zeros(M,N);
for n=1:nsum;
   Z=Z+(b(n)*sin(n*pi*x./l))'*exp(-((n*pi*al/l)^2)*t);
end;
Z=Z';
%%%%%%%% Matthewsin kokoelmasta: %%%%%%%%%%%%%%%%%%%%%%%%%
%
% Lämpöyhtälö
% Finite difference method 
% Sijoita tiedostoon yforwdif.m
%
function [U,x,t] = yforwdif(f,g1,g2,a,b,c,n,m)
% [U,x,t] = yforwdif(f,g1,g2,a,b,c,n,m)
% Finite difference solution to the heat equation.
% f  is a function, input.
% g1 is a function, input.
% g2 is a function, input.
% a  is the width of [0 a], input.
% b  is the width of [0 b], input.
% c  is the constant in the heat equation, input.
% n  is the number of grid points over [0 a], input.
% m  is the number of grid points over [0 b], input.
% U  is the solution matrix, output.
x=linspace(0,a,n);t=linspace(0,b,m);  % HA:n lis{ys
h = a/(n-1);
k = b/(m-1);
r = c^2*k/h^2;
s = 1 - 2*r;
U = zeros(n,m);
for j = 1:m,
   U(1,j) = feval(g1,k*(j-1));
   U(n,j) = feval(g2,k*(j-1));
end
for i = 2:(n-1),
   U(i,1) = feval(f,h*(i-1));
end
for j = 2:m,
  for i = 2:(n-1),
     U(i,j) = s*U(i,j-1) + r*(U(i-1,j-1) + U(i+1,j-1));
  end
end
U=U';
%%%%%%%%%%%% END OF yforwdif %%%%%%%%%%%%%%%%%%%%%%%%
%
%  Lämpöyhtälö: Crank-Nicholson  (implisiittinen differenssimenet.)
%  ycrnich  Matthewsin kokoelmasta
%  Sijoita tiedostoon ycrnich.m 
%  Tarvitsee myös alla olevan funktion trisys
%
function [U,x,t] = ycrnich(f,g1,g2,a,b,c,n,m)
% [U,x,t] = ycrnich(f,g1,g2,a,b,c,n,m)
% Crank-Nicholson solution to the heat equation.
% This subroutine uses trisys.m to solve a tri-diagonal system.
% f  is a function, input.
% g1 is a function, input.
% g2 is a function, input.
% a  is the width of [0 a], input.
% b  is the width of [0 b], input.
% c  is the constant in the heat equation, input.
% n  is the number of grid points over [0 a], input.
% m  is the number of grid points over [0 b], input.
% U  is the solution matrix, output.
x=linspace(0,a,n);t=linspace(0,b,m);  % HA:n lis{ys
h = a/(n-1);
k = b/(m-1);
r = c^2*k/h^2;
s1 = 2 + 2/r;
s2 = 2/r - 2;
U = zeros(n,m);
for j=1:m,
   U(1,j) = feval(g1,k*(j-1));
   U(n,j) = feval(g2,k*(j-1));
end
for i=2:(n-1),
   U(i,1) = feval(f,h*(i-1));
end
Vd = s1*ones(1,n);
Vd(1) = 1;
Vd(n) = 1;
Va = -ones(1,n-1);
Va(n-1) = 0;
Vc = -ones(1,n-1);
Vc(1) = 0;
Vb(1) = feval(g1,k*0);
Vb(n) = feval(g2,k*0);
for j=2:m,
  for i=2:(n-1),
     Vb(i) = U(i-1,j-1) + U(i+1,j-1) + s2*U(i,j-1);
  end
   X = trisys(Va,Vd,Vc,Vb);
   U(1:n,j) = X';
end
%%%%%%%%%%%%%%%% END of yrnich %%%%%%%%%%%%%%%%%%
% Käyttöesimerkki:
%   Tarvitset reuna- ja alkuarvofunktiot. Ota tästä ja modifioi tarpeen
%   mukaan
%%% 
%%% Tiedostoon reuna.m 
function z=reuna(t)
z=0;       
%%%%%%%%%%%%%%%%%%%
%%  Tiedostoon alku.m 
function z=alku(x)
z=100;       
%%%%%%%%%%%%%%%%%%%
% Nyt voit kutsua vaikka:
[U,x,t] = ycrnich('alku','reuna','reuna',1,10,1,5,25);
surfc(t,x,U)
%
% Matthewsin kokoelmasta: trisys (tarvitaan ainakin ycrnich:ssä)
% Sijoita tiedostoon trisys.m
%
function X = trisys(A,D,C,B)
% X = trisys(A,D,C,B)
% Solution of a triangular linear system AX = B.
% It is assumed that D and B have dimension n,
% and that A and C have dimension n-1;
% A sub diagonal, input.
% D diagonal vector, input.
% C super diagonal, input.
% B right hand side vector, input.
% X  solution vector, output.
n = length(B);
for k = 2:n,
  mult = A(k-1)/D(k-1);
  D(k) = D(k) - mult*C(k-1);
  B(k) = B(k) - mult*B(k-1);
end
X(n) = B(n)/D(n);
for k = (n-1):-1:1,
  X(k) = (B(k) - C(k)*X(k+1))/D(k);
end
%%%%%%%%%%%%%%% END of trisys %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Laplacen yhtälö
%  ydirich  Matthewsin kokoelmasta
%  Sijoita tiedostoon ydirich.m 
%
%%%%%%%%%%%% BEGIN ydirich %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U,x,t] = ydirich(f1,f2,f3,f4,a,b,h,tol,max1)
% [U,x,t] = ydirich(f1,f2,f3,f4,a,b,h,tol,max1)
% Dirichlet solution to Laplace's equation.
% f1  is a boundary function, input.
% f2  is a boundary function, input.
% f3  is a boundary function, input.
% f4  is a boundary unction, input.
% a   is the width of [0 a], input.
% b   is the width of [0 b], input.
% h   is the step size, input.
% tol is the tolerance, input.
% max1 is the maximum number of iterations, input.
% U  is the solution matrix, output.
n = fix(a/h)+1;
m = fix(b/h)+1;
x=linspace(0,a,n);t=linspace(0,b,m);  % HA:n lis{ys
ave = (a*(feval(f1,0)+feval(f2,0)) ...
    +  b*(feval(f3,0)+feval(f4,0)))/(2*a+2*b);
U = ave*ones(n,m);
for j=1:m,
  U(1,j) = feval(f3,h*(j-1));
  U(n,j) = feval(f4,h*(j-1));
end
for i=1:n,
  U(i,1) = feval(f1,h*(i-1));
  U(i,m) = feval(f2,h*(i-1));
end
U(1,1) = (U(1,2) + U(2,1))/2;
U(1,m) = (U(1,m-1) + U(2,m))/2;
U(n,1) = (U(n-1,1) + U(n,2))/2;
U(n,m) = (U(n-1,m) + U(n,m-1))/2;
w = 4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2));
err = 1;
cnt = 0;
while ((err>tol)&(cnt<=max1))
  err = 0;
  for j=2:(m-1),
    for i=2:(n-1),
      relx = w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+ U(i-1,j)-4*U(i,j))/4;
      U(i,j) = U(i,j) + relx;
      if (err<=abs(relx)), err=abs(relx); end
    end
  end
  cnt = cnt+1;
end
%%%%%%%%%%%%%%%% END of ydirich %%%%%%%%%%%%%%%%%%%%%%%%%%%