% matlab/osdynum.m  Luento 27.11.03
% Crank-Nicholson'n menetelmä lämpöyhtälölle
% Implisiittinen menetelmä, absoluuttisesti stabiili
% Esim. 1 KRE ss. 976-978
%
clear; close all
x=0:0.2:1
u0=sin(pi*x)       % Alkuarvot vektoriin u0
u0=u0';
n=5;
A=4*eye(n-1)+diag(-1*ones(n-2,1),-1)+diag(-1*ones(n-2,1),1)

%     1. aika-askel j:0->1

b=[u0(3);u0(2)+u0(4);u0(3)+u0(5);u0(4)]
u1=A\b;
u1=[0;u1;0];     % 1. aikataso vektorissa u1  t=k=0.04
U=[u0 u1]'

%%    2. aika-askel j:1->2

b=[u1(3);u1(2)+u1(4);u1(3)+u1(5);u1(4)]
u2=A\b;
u2=[0;u2;0];     % 2. aikataso vektorissa u2  t=2*k=0.08
U=[u0 u1 u2]'

%%%    3. aika-askel j:2->3

b=[u2(3);u2(2)+u2(4);u2(3)+u2(5);u2(4)]
u3=A\b;
u3=[0;u3;0];     % 3. aikataso vektorissa u3  t=3*k=0.12
U=[u0 u1 u2 u3]'

t=0:0.04:0.12
mesh(x,t,U)
figure
plot(x,U)
legend('t=0','t=0.04','t=0.08','t=0.12')
title('Crank-Nicholsonilla lasketun ratkaisun aikasiivuja')

% Verrataan tarkkaan ratkaisuun u(x,t)= sin(pi*x)*exp(-pi^2*t)

[X,T]=meshgrid(x,t)

Ut=sin(pi*X).*exp(-pi^2*T)     % Tarkka ratkaisu hilapisteistössä.

U,Ut                      % Taulukot allekkain tulostettuna

U-Ut                        % Ratkaisujen erotustaulukko

% tai pystyvektori*vaakavektori (oikea matriisikertolasku)
% exp(-pi^2*t)'*sin(pi*x)

mesh(x,t,Ut)

plot(x,Ut)
legend('t=0','t=0.02','t=0.04','t=0.06')
title('Tarkan ratkaisun aikasiivuja')

hold off
mesh(x,t,Ut); hold on;shg
surf(x,t,U-0.1); shg
