微分方程数值解2

Ps:题目取自李荣华第四版
交错方向迭代求解带边值的微分方程
一、实验目的及要求
Delta u= 2*pi^2*exp(pi*(x+y))*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y))
(x,y)属于G：（0,1）*（0,1）
U=0,(x,y)属于G的边界
其精确解为u(x,y)= exp(pi*(x+y))*sin(pi*x)*sin(pi*y)
取步长h=k=1/64,1/128，用五点差分格式和交替方向法求解。

二、代码和结果分析
（1）五点差分法
%p为精确解，e为误差，u为函数值，x，y为网格点
%h为步长，n*n为网格点数，kmax为最大迭代次数，ep为可控误差
clear
clc
syms temp;
n=64;
ep=1e-6;
kmax=1000;
u=zeros(n+1,n+1);
h=1/n;
x=0:h:1;
y=0:h:1;
%³õÖµ
for i=1:n+1
u(i,1)=0;
u(i,n+1)=0;
u(1,i)=0;
u(n+1,i)=0;
end
for i=1:n
for j=1:n
f(i,j)=2*pi^2*exp(pi*(x(i)+y(j)))*(sin(pi*x(i))*cos(pi*y(j))+cos(pi*x (i))*sin(pi*y(j)));
end
end
t=zeros(n-1,n-1);
for k=1:kmax
for i=2:n
for j=2:n
temp=-h^2*f(i,j)/4+(u(i,j+1)+u(i,j-1)+u(i+1,j)+u(i-1,j))/4; t(i,j)=(temp-u(i,j))^2;
u(i,j)=temp;
end
end
t(i,j)=sqrt(t(i,j));
if k>kmax
break;
end
if max(max(t))<ep
break;
end
end
for i=1:n+1
for j=1:n+1
p(i,j)=exp(pi*(x(i)+y(j)))*sin(pi*x(i))*sin(pi*y(j));
e(i,j)=abs(u(i,j)-p(i,j));
end
end
u
p
%五点差分求解图
figure;mesh(x,y,u)
%原图
figure;mesh(x,y,p)
%误差图
figure;mesh(x,y,e)
分析：从误差结果来看最大误差不超过3，差分求解结果比较接近真实值，进一步发现误差成正态分布，符合初边值问题的微分方程求解。

(2)交替方向迭代clc;
clear;
%N： X,Y的网格点
N=64+1;
% 进行网格剖分
n=N-1;
h=1/n;
%得到网格坐标
x=0:h:1;
y=0:h:1;
%初始化网格函数 u1 u2
u1=zeros(N,N);
u2=zeros(N,N);
for i=1:n
for j=1:n
f(i,j)=-2*pi^2*exp(pi*(x(i)+y(j)))*(sin(pi*x(i))*cos(pi*y(j))+cos(pi* x(i))*sin(pi*y(j)));
end
end
kmax=1000;
t=h^2/2;
for kk=1:kmax
for k=2:n
A=sparse(N,N);b=zeros(N,1);
for j=2:n
A(j,j-1)=-1/2;
A(j,j)=2;
A(j,j+1)=-1/2;
b(j)=1/2*u1(j,k-1)+1/2*u1(j,k+1)+t*f(j,k);
end
A(1,1)=1;
A(N,N)=1;
b(1)=0;
b(N)=0;
ut=A\b;
u2(:,k)=ut;
end
u1=u2;
for j=2:n
A=sparse(N,N);b=zeros(N,1);
for k=2:n
b(k)=1/2*u1(j-1,k)+1/2*u1(j+1,k)+t*f(j,k);
A(k,k-1)=-1/2;
A(k,k)=2;
A(k,k+1)=-1/2;
end
A(1,1)=1;
A(N,N)=1;
b(1)=0;
b(N)=0;
ut=A\b;
u2(j,:)=ut;
end
u1=u2;
end
figure;
mesh(1:N,1:N,u1)
e=zeros(N,N);
for i=1:N
for j=1:N
p(i,j)=exp(pi*(x(i)+y(j)))*sin(pi*x(i))*sin(pi*y(j));
e(i,j)=abs(u1(i,j)-p(i,j));
end
end
figure;mesh(1:N,1:N,e)
分析：交替方向迭代的误差不超过0.2，同样成正态分布。

三、总结
对比两种方法，在迭代次数相同的情况下，交替方向迭代的误差小于五点差分法。

在收敛速度上，通过比较等距离次数的迭代看误差减小幅度，也是交替方向迭代快点。