Metoda Crank - Nicolsona
: 14 sie 2022, o 11:51
Metoda Crank-Nicolsona
Metoda Crank-Nicolsona różni się od metody poprzedniej równaniem różnicowym, w którym pochodną cząstkową \(\displaystyle{ u_{xx},}\) zastąpiono średnią arytmetyczną jej różnic centralnych.
\(\displaystyle{ \frac{u_{i}^{j+1} - u_{i}^{j}}{k} = \frac{\alpha}{2}\left[ \frac{u_{i+1}^{j}-2u^{j}_{i}+u_{i-1}^{j}}{h^2} +\frac{u_{i+1}^{j+1}-2u^{j+1}_{i}+u_{i-1}^{j+1}}{h^2}\right]. }\)
Ta zmiana powoduje osiągnięcie większej dokładności wyników i stabilności numerycznej algorytmu.
Podstawiamy \(\displaystyle{ \lambda = \alpha \frac{k}{h^2}, }\) otrzymując schemat różnicowy:
\(\displaystyle{ -\lambda u_{i-1}^{j+1} +2(1+\lambda)u_{i}^{j+1}= \lambda u_{i-1}^{j}+2(1-\lambda)u_{i}^{j}+\lambda u_{i+1}^{j}, }\)
\(\displaystyle{ i=1,2,...,n-1, \ \ j = 0,1, \ \ ... m-1.}\)
Można pokazać, że metoda jest rzędu dokładności \(\displaystyle{ O(k^2 + h^2).}\)
Program w MATLAB
Rozwiążemy to samo równanie
\(\displaystyle{ u_{xx} = u_{t}, \ \ 0 < x < 1, \ \ , t >0 , }\)
z warunkami brzegowymi:
\(\displaystyle{ u(x,0) = \sin(\pi x) , \ \ 0\leq x \leq 1, \ \ u(0,t) = u(1,t) =0, \ \ t\geq 0 }\)
dla kroków dyskretyzacji różnej długości.
function f=f1(x,t);
f= sin(pi*x);
[/code]
Przykład 1
Dodano po 18 minutach 43 sekundach:
Przykład 2
Metoda Crank-Nicolsona różni się od metody poprzedniej równaniem różnicowym, w którym pochodną cząstkową \(\displaystyle{ u_{xx},}\) zastąpiono średnią arytmetyczną jej różnic centralnych.
\(\displaystyle{ \frac{u_{i}^{j+1} - u_{i}^{j}}{k} = \frac{\alpha}{2}\left[ \frac{u_{i+1}^{j}-2u^{j}_{i}+u_{i-1}^{j}}{h^2} +\frac{u_{i+1}^{j+1}-2u^{j+1}_{i}+u_{i-1}^{j+1}}{h^2}\right]. }\)
Ta zmiana powoduje osiągnięcie większej dokładności wyników i stabilności numerycznej algorytmu.
Podstawiamy \(\displaystyle{ \lambda = \alpha \frac{k}{h^2}, }\) otrzymując schemat różnicowy:
\(\displaystyle{ -\lambda u_{i-1}^{j+1} +2(1+\lambda)u_{i}^{j+1}= \lambda u_{i-1}^{j}+2(1-\lambda)u_{i}^{j}+\lambda u_{i+1}^{j}, }\)
\(\displaystyle{ i=1,2,...,n-1, \ \ j = 0,1, \ \ ... m-1.}\)
Można pokazać, że metoda jest rzędu dokładności \(\displaystyle{ O(k^2 + h^2).}\)
Program w MATLAB
Rozwiążemy to samo równanie
\(\displaystyle{ u_{xx} = u_{t}, \ \ 0 < x < 1, \ \ , t >0 , }\)
z warunkami brzegowymi:
\(\displaystyle{ u(x,0) = \sin(\pi x) , \ \ 0\leq x \leq 1, \ \ u(0,t) = u(1,t) =0, \ \ t\geq 0 }\)
dla kroków dyskretyzacji różnej długości.
Kod: Zaznacz cały
function crank_nicolson(f,c1,c2,L,T,h,k,alpha)
% Rozwiązanie równania ciepła z warunkami: u(x,0)=f(x)i u(0,t)=c1,
% i u(L,t)=c2 metodą Crank-Nicolsona.
n=L/h; m=T/k; lambda=alpha^2*k/(h^2)
z=0:h:L;
disp('_____________________________________________')
fprintf(' t x = ')
fprintf('%4.2f ',z)
fprintf('\n')
disp('______________________________________________')
fprintf('% 4.2f ',0)
for i=1:n+1
u(i)=feval(f,(i-1)*h);
fprintf('%10.6f ',u(i))
end
fprintf('\n')
for i=2:n
if (i~=n)
a(i)=-lambda;
end
b(i)=2*(1+lambda);
if (i~=n)
c(i)=-lambda;
end
end
bb=b;
for j=1:m
t=j*k;
fprintf('% 4.2f ',t)
for i=2:n
d(i)=lambda*u(i-1)+2*(1-lambda)*u(i)+lambda*u(i+1);
end
y(n+1)=c2;
y(1)=c1;
for i=3:n
ymult=a(i-1)/bb(i-1);
bb(i)=bb(i)-ymult*c(i-1);
d(i)=d(i)-ymult*d(i-1);
end
y(n)=d(n)/bb(n);
for i=n-1:-1:2
y(i)=(d(i)-c(i)*y(i+1))/bb(i);
end
for i=1:n+1
fprintf('%10.6f ',y(i))
end
fprintf('\n')
u=y;
bb=b;
end
f= sin(pi*x);
[/code]
Przykład 1
Kod: Zaznacz cały
>> crank_nicolson('f1',0,0,1,0.5,0.2,0.05,1)
lambda =
1.2500
________________________________________________________________________
t\ x = 0.00 0.20 0.40 0.60 0.80 1.00
________________________________________________________________________
0.00 0.000000 0.587785 0.951057 0.951057 0.587785 0.000000
0.05 0.000000 0.361228 0.584480 0.584480 0.361228 0.000000
0.10 0.000000 0.221996 0.359197 0.359197 0.221996 0.000000
0.15 0.000000 0.136430 0.220748 0.220748 0.136430 0.000000
0.20 0.000000 0.083844 0.135662 0.135662 0.083844 0.000000
0.25 0.000000 0.051527 0.083372 0.083372 0.051527 0.000000
0.30 0.000000 0.031666 0.051237 0.051237 0.031666 0.000000
0.35 0.000000 0.019461 0.031488 0.031488 0.019461 0.000000
0.40 0.000000 0.011960 0.019351 0.019351 0.011960 0.000000
0.45 0.000000 0.007350 0.011893 0.011893 0.007350 0.000000
0.50 0.000000 0.004517 0.007309 0.007309 0.004517 0.000000
Przykład 2
Kod: Zaznacz cały
crank_nicolson('f1',0,0,1,.025,.1,.0025,1)
lambda =
0.2500
______________________________________________________________________________________________________________________________
t\ x = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
______________________________________________________________________________________________________________________________
0.00 0.000000 0.309017 0.587785 0.809017 0.951057 1.000000 0.951057 0.809017 0.587785 0.309017 0.000000
0.01 0.000000 0.301546 0.573575 0.789458 0.928064 0.975824 0.928064 0.789458 0.573575 0.301546 0.000000
0.02 0.000000 0.294256 0.559708 0.770372 0.905627 0.952233 0.905627 0.770372 0.559708 0.294256 0.000000
0.03 0.000000 0.287142 0.546177 0.751748 0.883733 0.929212 0.883733 0.751748 0.546177 0.287142 0.000000
0.04 0.000000 0.280200 0.532972 0.733574 0.862368 0.906747 0.862368 0.733574 0.532972 0.280200 0.000000
0.05 0.000000 0.273426 0.520087 0.715839 0.841519 0.884826 0.841519 0.715839 0.520087 0.273426 0.000000
0.06 0.000000 0.266816 0.507514 0.698533 0.821175 0.863434 0.821175 0.698533 0.507514 0.266816 0.000000
0.07 0.000000 0.260365 0.495244 0.681645 0.801322 0.842560 0.801322 0.681645 0.495244 0.260365 0.000000
0.08 0.000000 0.254071 0.483271 0.665166 0.781949 0.822190 0.781949 0.665166 0.483271 0.254071 0.000000
0.09 0.000000 0.247928 0.471588 0.649085 0.763045 0.802313 0.763045 0.649085 0.471588 0.247928 0.000000
0.10 0.000000 0.241934 0.460187 0.633392 0.744598 0.782916 0.744598 0.633392 0.460187 0.241934 0.000000