% Viimeinen luento 9.12.05

% Luentoesimerkki: Neliö [0,4]x[0,4], AE ylhäällä: u(x,4)=f(x)=1/4*x*(4-x)
clear; close all  
fplot('1/4*x*(4-x)',[0 4])
title('Lämpötilajakauma yläreunalla')

c=inline('-16*(-1+(-1)^n)/n^3/pi^3/sinh(n*pi)','n')
a=4;b=4;
M=20;N=20;nb=20;
x=linspace(0,a,M); y=linspace(0,b,N);
Z=zeros(M,N);
for n=1:nb;
   Z=Z+c(n)*sin(n*pi*x/a)'*sinh(n*pi*y/a);
end;
Z=Z';

figure
surfc(x,y,Z)
colorbar
title('Analyyttinen (Fourir-sarja-)ratkaisu')
% Ratkaistaan sama differenssimenetelmällä.
% Aluksi rakennetaan käsin pieni: 4=2x2 sisäsolmua, eli väli [0,4] jaetaan 3:een 
% osaan
 
 
f=inline('1/4*x.*(4-x)','x')
%fplot(f,[0 4])
%title('Lämpötilajakauma yläreunalla')

h=4/3; x=0:h:4;   % x- ja y-diskretointi.
% Reuna-arvovektorit:
E=[0 0 0 0];     % Etelä
L=E';Ita=L;      % L -- länsi, Ita on sama, laitetaan huvin vuoksi pystyyn.
P=f(x)           % Pohjoinen: lasketaan f:n arvot solmupisteissä.

% Rakennetaan Laplace-diskretointimatriisi A:
B=[-4 1;1 -4]
I=eye(2,2)
A=[B I;I B]   % Lohkorakenne.
b=-[E(2)+L(2);E(3)+Ita(2);L(3)+P(2);Ita(3)+P(3)]
u=A\b         % Ratkaisuvektori
U=reshape(u,2,2)
U=flipud(U');   % Muotoillaan sisäsolmumuotoon.

% Kehystetään reuna-arvoilla:
W=[L(2:3) U Ita(2:3)];      
W=[P;W;E];

%W(1,1)=(W(1,2)+W(2,1))/2;   % Ei näitä, lapsolve1 tekee, mutta tässä vain haittaa.
%W(1,4)=(W(1,3)+W(4,2))/2

figure
surfc(x,x,flipud(W)); colorbar; title('Differenssimenetelmä, 2 x 2 sisäsolmua');shg
%figure
%contour(x,x,flipud(W),20)

% Kokeillaan aiemmilla kurseilla tehtyä valmista lapsolve1-ratkaisjaamme:

Z=lapsolve1(E(2:3),Ita(2:3)',P(2:3),L(2:3)')
Z-W   % Sisäsolmujen osalta samat, reunat käsitellään lapsolvessa keskiarvoistaen.
      % Tämän tehtävän kannalta oikeampi on edellä oleva tapa. (Kun h:ta pienennetään
      % ero mitätöityy.

      
% 4 sisäsolmua:

h=4/5; x=0:h:4;
ala=[0 0 0 0]; vasen=ala; oikea=ala;
yla=f(x); yla=yla(2:end-1)
U4=lapsolve1(ala,oikea,yla,vasen)

[U4,A]=lapsolve1(ala,oikea,yla,vasen)
figure
surfc(x,x,U4); colorbar; title('Differenssimenetelmä, 4 x 4 sisäsolmua');shg

% Laplace-diskretointimatriisi on nyt tämännäköinen.
full(A)

% Tästä voisi jatkaa tutkimalla, miten samoissa hilapisteissä laskettujen arvojen
% virhe käyttäytyy, kun jakoa hienonnetaan. (Kun solmumäärä koordinaattisuunnassa
% kaksinkertaistuu (tasossa nelinkertaistuu), niin miten virhe?. 
% No, taitaa tulla joulu vastaan! 

figure
spy(A);title('A-matriisin (16 x 16) rakenne spy(A)')


