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

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

摘要生成于 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;
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值