http://www.math.hut.fi/teaching/v/matlab/opas/esim/    Päivitetty 19.10.2001

Laplacen operaattorin diskretointi (5 pisteen klassinen)

Laplacen / Poissonin yhtälön delta u = 0 / delta u = f diskretointi 5 pisteen kaavalla johtaa differenssikaavaan: 
ui-1,j+ui,j-1-4 ui,j+ui,j+1+ui+1,j= h2 fi,j
Laplacen yhtälön tapauksessa siis oikea puoli on = 0.

y1y2 y3 y4
v4 u4,1 u4,2 u4,3 u4,4 o4
v3 u3,1 u3,2 u3,3 u3,4 o3
*** Jatka tästä taulukon editointia ***

v2 u2,1 u2,2 u2,3 u2,4 o2 v1 | u1,1 u1,2 u1,3 u1,4 o1 a1 a2 a3 a4 


      50    20    10      0
     |-----|-----|-----|------|-----|
     |                              |
  80    | u[4,1] u[4,2] u[4,3] u[4,4]  | 20
     |                              |
     |                              |
  60    | u[3,1] u[3,2] u[3,3] u[3,4]  | 20
     |                              |
     |                              |
  40    | u[2,1] u[2,2] u[2,3] u[2,4]  | 20
     |                              |
     |                              |
  20    | u[1,1] u[1,2] u[1,3] u[1,4]  | 20
     |                              |
     |-----|-----|-----|------|-----|
     0     0     0      0
Kun edetään vasemmalta oikealle ja alhaalta ylös, eli indeksi lasketaan kaavalla (i,j) -> n(i-1)+j , saadaan yhtälösysteemin matriisi: 

       [-4 , 1 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0]
       [                                                              ]
       [1 , -4 , 1 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0]
       [                                                              ]
       [0 , 1 , -4 , 1 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0]
       [                                                              ]
       [0 , 0 , 1 , -4 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0]
       [                                                              ]
       [1 , 0 , 0 , 0 , -4 , 1 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0]
       [                                                              ]
       [0 , 1 , 0 , 0 , 1 , -4 , 1 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0]
       [                                                              ]
       [0 , 0 , 1 , 0 , 0 , 1 , -4 , 1 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0]
       [                                                              ]
       [0 , 0 , 0 , 1 , 0 , 0 , 1 , -4 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0]
       [                                                              ]
       [0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , -4 , 1 , 0 , 0 , 1 , 0 , 0 , 0]
       [                                                              ]
       [0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 1 , -4 , 1 , 0 , 0 , 1 , 0 , 0]
       [                                                              ]
       [0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 1 , -4 , 1 , 0 , 0 , 1 , 0]
       [                                                              ]
       [0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 1 , -4 , 0 , 0 , 0 , 1]
       [                                                              ]
       [0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , -4 , 1 , 0 , 0]
       [                                                              ]
       [0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 1 , -4 , 1 , 0]
       [                                                              ]
       [0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 1 , -4 , 1]
       [                                                              ]
       [0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 1 , -4]
Yksityiskohtia seuraavissa: Skannatuja kalvoja (epäesteettisiä)
Ratkaisuohjeita kurssisivulta    Samoin 

    function A=lapmatr(n)
    % LAPMATR       A=lapmatr(n)
    %               This function forms matrix A for the numerical         
    %               approximation of the solution of Laplace equation. 
    %               The desired dimension n i.e the number of grid points 
    %               is given as a parameter.
    %
    N=n^2;
    maind=ones(1,N);
    fardiag=ones(1,N-n);
    j=1:N-1;  
    neardiag=rem(j,n)~=0;
    A=4*diag(maind)-diag(neardiag,1)-diag(neardiag,-1)- ...
    diag(fardiag,n)-diag(fardiag,-n);

function u=lapsolve(bot,right,top,left) 
    %
    % LAPSOLVE      u=lapsolve(bot,right,top,left)
    %               Solve Laplace equation on square. Parameters bot,right,top,
    %               left are vectors of boundary values. The vertical boundaries 
    %               are given from bottom to top and the horizontal boundaries
    %               from left to right. 
    n=length(bot);
    A=lapmatr(n);
    B=[top;zeros((n-2),n);bot]+[fliplr(left)',zeros(n,(n-2)),fliplr(right)'];
    B2=B';
    Tvek=A\B2(:);
    T=reshape(Tvek,n,n)';
    NW=0.5*(left(n)+top(1));NE=0.5*(top(n)+right(n));
    SE=0.5*(bot(n)+right(1));SW=0.5*(bot(1)+left(1));
    T2=zeros(n+2,n+2);
    T2(:,1)=[0; fliplr(left)'; 0];
    T2(:,n+2)=[0;  fliplr(right)'; 0];
    T2(1,:)=[NW top NE];
    T2(n+2,:)=[SW bot SE];
    T2([2:n+1],[2:n+1])=T;
    u=T2;

function [bot,right,top,left]=bound(n)
    %%
    %%BOUND [bot,right,top,left]=bound(n)
    %% 
    %%
    p3=polyfit([0,0.25,0.5 0.75],[50 20 10 0],3);
    p1=polyfit([0 0.25 0.5 0.75],[20 40 60 80],1);
    x=0:(0.75/n):0.75; y=x;
    top=polyval(p3,x);bot=zeros(1,length(top));
    right=polyval(p1,y); left=right;


    >> [bot,right,top,left]=bound(30);
    >> T=lapsolve(bot,right,top,left);
    >> subplot(2,2,1);mesh(u);view(-37.5,20);
    >> subplot(2,2,2);mesh(u);view(-37.5+90,20);
    >> subplot(2,2,3);mesh(u);view(-37.5+180,20);
    >> subplot(2,2,4);mesh(u);view(-37.5+270,20);