% matlab/lampo/CRNich.m  Luento 5.12.08
% 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
n=length(x)  % Osavälien lkm.
u0=sin(pi*x)'       % Alkuarvot vektoriin u0
u0s=u0(2:n-1)        % Sisäsolmuarvot
%A=4*eye(n-1)+diag(-1*ones(n-2,1),-1)+diag(-1*ones(n-2,1),1)
% Tehokkampaa ja helpompi kirjoittaa spdiags:lla:
e=ones(n-2,1);
A=spdiags([-e 4*e -e],-1:1,n-2,n-2);
KatsotaanAtaysi=full(A)
B=spdiags([e e],[-1 1],4,4);
KatsotaanBtaysi=full(B)

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

%  Tässä olis vähemmän elegantti, vaikeammin yleistyvä tapa:
% u0=sin(pi*x)'; b=[u0(3);u0(2)+u0(4);u0(3)+u0(5);u0(4)];
%  Täsä se elegantti (vrt. luennot)
b=B*u0s
u1s=A\b;
u1=[0;u1s;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)]
b=B*u1s
u2s=A\b;
u2=[0;u2s;0]     % 2. aikataso vektorissa u2  t=2*k=0.08
U=[u0 u1 u2]'


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

b=B*u2s
u3s=A\b;
u3=[0;u3s;0];     % 3. aikataso vektorissa u3  t=3*k=0.12
U=[u0 u1 u2 u3]'

t=0:0.04:0.12
surf(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)
% Tehdään vasta ensi viikolla.
