Contents
% 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]'
x =
0 0.2000 0.4000 0.6000 0.8000 1.0000
n =
6
u0 =
0
0.5878
0.9511
0.9511
0.5878
0.0000
u0s =
0.5878
0.9511
0.9511
0.5878
KatsotaanAtaysi =
4 -1 0 0
-1 4 -1 0
0 -1 4 -1
0 0 -1 4
KatsotaanBtaysi =
0 1 0 0
1 0 1 0
0 1 0 1
0 0 1 0
b =
0.9511
1.5388
1.5388
0.9511
u1 =
0
0.3993
0.6460
0.6460
0.3993
0
U =
0 0.5878 0.9511 0.9511 0.5878 0.0000
0 0.3993 0.6460 0.6460 0.3993 0
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]'
b =
0.6460
1.0453
1.0453
0.6460
u2 =
0
0.2712
0.4388
0.4388
0.2712
0
U =
0 0.5878 0.9511 0.9511 0.5878 0.0000
0 0.3993 0.6460 0.6460 0.3993 0
0 0.2712 0.4388 0.4388 0.2712 0
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.
b =
0.4388
0.7101
0.7101
0.4388
U =
0 0.5878 0.9511 0.9511 0.5878 0.0000
0 0.3993 0.6460 0.6460 0.3993 0
0 0.2712 0.4388 0.4388 0.2712 0
0 0.1842 0.2981 0.2981 0.1842 0
t =
0 0.0400 0.0800 0.1200
