Matlab:Crank Nicolson方法求解线性抛物方程

最新推荐文章于 2024-04-10 16:44:52 发布

转载最新推荐文章于 2024-04-10 16:44:52 发布 · 2.1k 阅读

27 ·

CC 4.0 BY-SA版权

原文链接：https://www.cnblogs.com/xtu-hudongdong/p/6560734.html

文章标签：

#matlab

数值分析从入门到入土专栏收录该内容

3 篇文章

订阅专栏

该博客主要探讨了克兰克-尼科尔森方法在数值求解偏微分方程中的应用。通过设置不同的网格步数，计算并比较了数值解与精确解的误差，展示了误差随网格步数增加而减小的趋势。同时，通过三维图像分别呈现了精确解、数值解和误差分布，进一步分析了误差特性。最后，计算了误差阶数，验证了方法的稳定性。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

在这里插入图片描述

tic;
clear
clc
M=[10,20,40,80,160,320,640];%x的步数
K=M; %时间t的步数
for p=1:length(M)
hx=1/M(p);
ht=1/K(p);
r=ht/hx^2; %网格比
x=0:hx:1;
t=0:ht:1;
numerical=zeros(M(p)+1,K(p)+1);
numerical(:,1)=exp(x); %初始值
numerical(1,:)=exp(t); %边值
numerical(M(p)+1,:)=exp(t+1); %边值
a=-r/2*ones(M(p)-2,1);b=(1+r)*ones(M(p)-1,1);c=-r/2*ones(M(p)-2,1);
fun1=inline('exp(x+t)','x','t');
for i=1:length(x)
    for j=1:length(t)
Accurate(i,j)=fun1(x(i),t(j));
    end
end
d=r/2*ones(M(p)-2,1);e=(1-r)*ones(M(p)-1,1);f=r/2*ones(M(p)-2,1);
B=diag(d,-1)+diag(e,0)+diag(f,1);
fun2=inline('0','x','t');
for i=1:M(p)-1
for k=1:K(p)
    f(i,k)=ht*fun2(x(i+1),t(k)+ht/2);
end
end
for k=1:K(p)
f(1,k)=r/2*(numerical(1,k+1)+numerical(1,k));
f(M(p)-1,k)=r/2*(numerical(M(p)+1,k+1)+numerical(M(p)+1,k));
end
for k=1:K(p)
    right_vector=f(:,k)+B*numerical(2:M(p),k);
    numerical(2:M(p),k+1)=chase(a,b,c,right_vector);
end
error=numerical(2:M(p),2:K(p))'-Accurate(2:M(p),2:K(p))';
error_inf(p)=max(max(error));
figure(p)
[X,Y]=meshgrid(x,t);
subplot(1,3,1)
mesh(X,Y,Accurate');
xlabel('x'),ylabel('t');zlabel('Accurate');
title('the image of Accurate result');
grid on
subplot(1,3,2)
mesh(X,Y,numerical');
xlabel('x'),ylabel('t');zlabel('numerical');
title('the image of numerical result');
grid on
subplot(1,3,3)
mesh(X,Y,numerical'-Accurate');
xlabel('x'),ylabel('t');zlabel('error');
title('the image of error result');
grid on
end
for k=2:length(M)
    H=error_inf(p-1)/error_inf(p);
E_inf(k-1)=log2(H);
end
figure(length(M)+1)
plot(1:length(M)-1,E_inf,'-r v');
ylabel('误差阶数');
title('Crank nicolson 误差阶数');
grid on
toc;