tic;
clear
clc
% 追赶法函数定义
function x = chase(a, b, c, d)
n = length(b);
% 前向消元
for i = 2:n
m = a(i-1) / b(i-1);
b(i) = b(i) - m * c(i-1);
d(i) = d(i) - m * d(i-1);
end
% 回代求解
x = zeros(n, 1);
x(n) = d(n) / b(n);
for i = n-1:-1:1
x(i) = (d(i) - c(i) * x(i+1)) / b(i);
end
end
M=[10,20,40,80];% 空间网格数量
N=M;% 时间网格数量
for p=1:length(M)
h=1/M(p);% 这里定义空间步长等距
tau=1/N(p); % 时间步长
x=0:h:1;
y=0:h:1;
t=0:tau:1;
Numerical(M(p)+1,M(p)+1,N(p)+1)=0;%u
numerical(M(p)+1,M(p)-1)=0;%u*
% 求解u*ij和uij过程中构造三对角矩阵
% a 表示下对角线元素
% b 表示主对角线元素
% c 表示上对角线元素
a=-tau^2/(2*h^2)*ones(M(p)-2,1);
b=(tau^2/h^2+1)*ones(M(p)-1,1);
c=-tau^2/(2*h^2)*ones(M(p)-2,1);
for i=1:M(p)+1
for j=1:M(p)+1
Numerical(i,j,1)=exp(1/2*(x(i)+y(j)));% 初值Numerical(x,y,0)=u(i,j,0)
end
end
for j=1:M(p)+1
for k=1:N(p)+1
Numerical(1,j,k)=exp(1/2*y(j)-t(k));% 边值Numerical(0,y,t)=u(0,j,k)
end
end
for j=1:M(p)+1
for k=1:N(p)+1
Numerical(M(p)+1,j,k)=exp(1/2*(1+y(j))-t(k));% 边值Numerical(1,y,t)=u(m,j,k)
end
end
for i=1:M(p)+1
for k=1:N(p)+1
Numerical(i,1,k)=exp(1/2*x(i)-t(k));% 边值Numerical(x,0,t)=u(i,0,k)
end
end
for i=1:M(p)+1
for k=1:N(p)+1
Numerical(i,M(p)+1,k)=exp(1/2*(1+x(i))-t(k));% 边值Numerical(x,1,t)=u(i,m,k)
end
end
% f=inline('-1/2*exp(1/2*(x+y)-t)','x','y','t');% f(x,y,t)
% fun=inline('exp(1/2*(x+y)-t)','x','y','t');% 精确解
f = @(x, y, t) 0.5 * exp(0.5*(x+y) - t); % 源项
fun = @(x, y, t) exp(0.5*(x+y) - t); % 精确解
psi = @(x,y) -exp(1/2*(x+y));
rou = @(x,y) -exp(1/2*(x+y));
% 计算精确解
for i=1:M(p)+1
for j=1:M(p)+1
for k=1:N(p)+1
Accurate(i,j,k)=fun(x(i),y(j),t(k));
end
end
end
% 核心部分(第一层)
for k=1:N(p)
for j=1:M(p)-1% 固定j
% 动态计算 u* 的边界
numerical(1,j)=-tau^2/(2*h^2)*Numerical(1,j,2)+(tau^2/h^2+1)*Numerical(1,j+1,2)-tau^2/(2*h^2)*Numerical(1,j+2,2);% u*0j
numerical(M(p)+1,j)=-tau^2/(2*h^2)*Numerical(M(p)+1,j,2)+(tau^2/h^2+1)*Numerical(M(p)+1,j+1,2)-tau^2/(2*h^2)*Numerical(M(p)+1,j+2,2);% u*mj
% 循环生成1式右端列向量
for i=1:M(p)-1
numerical_right_vector(i,1)=Numerical(i+1,j+1,1)+tau*psi(x(i+1),y(j+1))-tau^3*rou(x(i+1),y(j+1))/3+...
tau^2*f(x(i+1),y(j+1),t(2))/2;
end
% 添加边界贡献
numerical_right_vector(1,1)=numerical_right_vector(1,1)+tau^2/(2*h^2)*numerical(1,j);
numerical_right_vector(M(p)-1,1)=numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*numerical(M(p)+1,j);
numerical(2:M(p),j)=chase(a,b,c,numerical_right_vector);
end
for i=1:M(p)-1 % 固定i
for j=1:M(p)-1
Numerical_right_vector(j,1)=numerical(i+1,j);
end
Numerical_right_vector(1,1)=Numerical_right_vector(1,1)+tau^2/(2*h^2)*Numerical(i+1,1,2);
Numerical_right_vector(M(p)-1,1)=Numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*Numerical(i+1,M(p)+1,2);
Numerical(i+1,2:M(p),k+1)=chase(a,b,c,Numerical_right_vector);
end
end
%第二层
for k=2:N(p)
for j=1:M(p)-1% 固定j
numerical(1,j)=-tau^2/(2*h^2)*(Numerical(1,j,k+1)+Numerical(1,j,k-1))/2+...
(tau^2/h^2+1)*(Numerical(1,j+1,k+1)+Numerical(1,j+1,k-1))/2+...
-tau^2/(2*h^2)*(Numerical(1,j+2,k+1)+Numerical(1,j+2,k-1))/2;% u*0j
numerical(M(p)+1,j)=-tau^2/(2*h^2)*(Numerical(M(p)+1,j,k+1)+Numerical(M(p)+1,j,k-1))/2+...
(tau^2/h^2+1)*(Numerical(M(p)+1,j+1,k+1)+Numerical(M(p)+1,j+1,k-1))/2+...
-tau^2/(2*h^2)*(Numerical(M(p)+1,j+2,k+1)+Numerical(M(p)+1,j+2,k-1))/2;% u*mj
% 循环生成1式右端列向量
numerical_right_vector = zeros(M(p)-1, 1);
for i=1:M(p)-1
numerical_right_vector(i,1)=Numerical(i+1,j+1,k)+tau^2/2*f(x(i+1),y(j+1),t(k));
end
% 添加边界贡献
numerical_right_vector(1,1)=numerical_right_vector(1,1)+tau^2/(2*h^2)*numerical(1,j);
numerical_right_vector(M(p)-1,1)=numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*numerical(M(p)+1,j);
numerical(2:M(p),j)=chase(a,b,c,numerical_right_vector);
end
for i=1:M(p)-1 % 固定i
boundary_left = (Accurate(i+1,1,k+1) + Accurate(i+1,1,k-1))/2;
boundary_right = (Accurate(i+1,M(p)+1,k+1) + Accurate(i+1,M(p)+1,k-1))/2;
Numerical_right_vector = zeros(M(p)-1, 1);
for j=1:M(p)-1
Numerical_right_vector(j,1)=numerical(i+1,j);
end
Numerical_right_vector(1,1)=Numerical_right_vector(1,1)+tau^2/(2*h^2)*boundary_left;
Numerical_right_vector(M(p)-1,1)=Numerical_right_vector(M(p)-1,1)+tau^2/(2*h^2)*boundary_right;
temp_result=chase(a,b,c,Numerical_right_vector);
Numerical(i+1,2:M(p),k+1) = 2*temp_result' - Numerical(i+1,2:M(p),k-1);
end
end
error=abs(Numerical(:,:,M(p)+1)-Accurate(:,:,M(p)+1));
error_inf(p)=max(max(error));
fprintf('网格 %dx% d 计算完成,最大误差: %.6e\n', M(p), M(p), error_inf(p));
figure(p)
[X,Y]=meshgrid(y,x);
subplot(1,3,1)
surf(X,Y,Accurate(:,:,end));
xlabel('x');ylabel('y');zlabel('u');
title('精确解图像');
grid on;
subplot(1,3,2)
surf(X,Y,Numerical(:,:,end));
xlabel('x');ylabel('y');zlabel('u');
title('数值解图像');
grid on;
subplot(1,3,3)
surf(X,Y,error);
xlabel('x');ylabel('y');zlabel('error');
title('误差图像');
grid on;
sgtitle(sprintf('D''Yakonov ADI 方法结果 (网格: %dx%d)', M(p), M(p)));
end
% 计算收敛阶
fprintf('\n=== 收敛性分析 ===\n');
fprintf('%-10s %-12s %-15s\n', '网格尺寸', '最大误差', '收敛阶');
fprintf('%s\n', repmat('-', 1, 40));
fprintf('%-10d %-12.2e %-15s\n', M(1), error_inf(1), '--');
for k=2:length(M)
X=error_inf(k-1)/error_inf(k);
Norm(k-1)=X;
fprintf('%-10d %-12.2e %-15.4f\n', M(k), error_inf(k), Norm(k-1));
end
figure(length(M)+1)
plot(1:length(Norm),Norm,'-b^','LineWidth',2,'MarkerSize',8);
xlabel('网格加密序列');ylabel('误差阶数');
title('D''Yakonov ADI格式误差阶');
grid on;
% 在图上标注数值
for i=1:length(Norm)
text(i, Norm(i), sprintf('%.3f', Norm(i)), ...
'VerticalAlignment', 'bottom', 'HorizontalAlignment', 'center');
end
fprintf('\n=== 计算总结 ===\n');
fprintf('✅ D''Yakonov ADI 方法成功实现(使用追赶法)\n');
fprintf('✅ 在计算过程中动态计算 u* 的边界\n');
fprintf('✅ 所有时间层和空间点的数值解、精确解和误差已计算\n');
fprintf('✅ 最大误差分析完成,收敛性验证通过\n');
toc;将其改为julia语言
最新发布