Metoda Crank - Nicolsona

Równania różniczkowe i całkowe. Równania różnicowe. Transformata Laplace'a i Fouriera oraz ich zastosowanie w równaniach różniczkowych.
janusz47
Użytkownik
Użytkownik
Posty: 7910
Rejestracja: 18 mar 2009, o 16:24
Płeć: Mężczyzna
Podziękował: 30 razy
Pomógł: 1670 razy

Metoda Crank - Nicolsona

Post autor: janusz47 »

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.

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 
function f=f1(x,t);
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 
 
 
Dodano po 18 minutach 43 sekundach:
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 
 
ODPOWIEDZ