Matlab:Crank Nicolson方法求解线性抛物方程
1 tic; 2 clear 3 clc 4 M=[10,20,40,80,160,320,640];%x的步数 5 K=M; %时间t的步数 6 for p=1:length(M) 7 hx=1/M(p); 8 ht=1/K(p); 9 r=ht/hx^2; %网格比 10 x=0:hx:1; 11 t=0:ht:1; 12 numerical=zeros(M(p)+1,K(p)+1); 13 numerical(:,1)=exp(x); %初始值 14 numerical(1,:)=exp(t); %边值 15 numerical(M(p)+1,:)=exp(t+1); %边值 16 a=-r/2*ones(M(p)-2,1);b=(1+r)*ones(M(p)-1,1);c=-r/2*ones(M(p)-2,1); 17 fun1=inline('exp(x+t)','x','t'); 18 for i=1:length(x) 19 for j=1:length(t) 20 Accurate(i,j)=fun1(x(i),t(j)); 21 end 22 end 23 d=r/2*ones(M(p)-2,1);e=(1-r)*ones(M(p)-1,1);f=r/2*ones(M(p)-2,1); 24 B=diag(d,-1)+diag(e,0)+diag(f,1); 25 fun2=inline('0','x','t'); 26 for i=1:M(p)-1 27 for k=1:K(p) 28 f(i,k)=ht*fun2(x(i+1),t(k)+ht/2); 29 end 30 end 31 for k=1:K(p) 32 f(1,k)=r/2*(numerical(1,k+1)+numerical(1,k)); 33 f(M(p)-1,k)=r/2*(numerical(M(p)+1,k+1)+numerical(M(p)+1,k)); 34 end 35 for k=1:K(p) 36 right_vector=f(:,k)+B*numerical(2:M(p),k); 37 numerical(2:M(p),k+1)=chase(a,b,c,right_vector); 38 end 39 error=numerical(2:M(p),2:K(p))'-Accurate(2:M(p),2:K(p))'; 40 error_inf(p)=max(max(error)); 41 figure(p) 42 [X,Y]=meshgrid(x,t); 43 subplot(1,3,1) 44 mesh(X,Y,Accurate'); 45 xlabel('x'),ylabel('t');zlabel('Accurate'); 46 title('the image of Accurate result'); 47 grid on 48 subplot(1,3,2) 49 mesh(X,Y,numerical'); 50 xlabel('x'),ylabel('t');zlabel('numerical'); 51 title('the image of numerical result'); 52 grid on 53 subplot(1,3,3) 54 mesh(X,Y,numerical'-Accurate'); 55 xlabel('x'),ylabel('t');zlabel('error'); 56 title('the image of error result'); 57 grid on 58 end 59 for k=2:length(M) 60 H=error_inf(p-1)/error_inf(p); 61 E_inf(k-1)=log2(H); 62 end 63 figure(length(M)+1) 64 plot(1:length(M)-1,E_inf,'-r v'); 65 ylabel('误差阶数'); 66 title('Crank nicolson 误差阶数'); 67 grid on 68 toc;
胡冬冬
