一维抛物线型方程数值解法(1)(附图及matlab程序)

一维抛物线偏微分方程数值解法（1）
解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法）
Ut-Uxx=0, 0<x<1,0<t<=1（Ut-aUxx=f(x,t),a>0)
U(x,0)=e^x, 0<=x<=1,
U(0,t)=e^t,U(1,t)=e^(1+t), 0<t<=1
精确解为：U(x,t)=e^(x+t);
下面给出两个matlab程序，实质一样（用的是向前欧拉格式）
第二个程序由之前解线性方程组的G-S迭代法得到，迭代次数k=2（固定）
function [p u e x t]=pwxywxq(h1,h2,m,n)
% 解抛物线型一维方程向前欧拉格式（Ut-aUxx=f(x,t),a>0)
%不用解线性方程组，由下一层(时间层）的值就直接得到上一层的值
%m,n为x,t方向的网格数，例如（2-0）/0.01=200;
%e为误差，p为精确解
u=zeros(n+1,m+1);
x=0+(0:m)*h1;
t=0+(0:n)*h2;
for(i=1:n+1)
u(i,1)=exp(t(i));
u(i,m+1)=exp(1+t(i));
end
for(i=1:m+1)
u(1,i)=exp(x(i));
end
for(i=1:n+1)
for(j=1:m+1)
f(i,j)=0;
end
end
r=h2/(h1*h1); %此处r=a*h2/(h1*h1）；a=1 要求r<=1/2差分格式才稳定for(i=1:n)
for(j=2:m)
u(i+1,j)=(1-2*r)*u(i,j)+r*(u(i,j-1)+u(i,j+1))+h2*f(i,j);
end
end
for(i=1:n+1)
for(j=1:m+1)
p(i,j)=exp(x(j)+t(i));
e(i,j)=abs(u(i,j)-p(i,j));
end
end
或者：
function [u e p x t k]=paowuxianyiweixq(h1,h2,m,n,kmax,ep)
% 解抛物线型一维方程向前欧拉格式（Ut-aUxx=f(x,t),a>0)
%kmax为最大迭代次数
%m,n为x,t方向的网格数，例如（2-0）/0.01=200;
%e为误差，p为精确解
syms temp;
u=zeros(n+1,m+1);
x=0+(0:m)*h1;
t=0+(0:n)*h2;
for(i=1:n+1)
u(i,1)=exp(t(i));
u(i,m+1)=exp(1+t(i));
end
for(i=1:m+1)
u(1,i)=exp(x(i));
end
for(i=1:n+1)
for(j=1:m+1)
f(i,j)=0;
end
end
a=zeros(n,m-1);
r=h2/(h1*h1);%此处r=a*h2/(h1*h1）；a=1 要求r<=1/2差分格式才稳定
for(k=1:kmax)
for(i=1:n)
for(j=2:m)
temp=(1-2*r)*u(i,j)+r*(u(i,j-1)+u(i,j+1))+h2*f(i,j); a(i+1,j)=(temp-u(i+1,j))*(temp-u(i+1,j));
u(i+1,j)=temp;
end
end
a(i,j)=sqrt(a(i,j));
if(k>kmax)
break;
end
if(max(max(a))<ep)
break;
end
end
for(i=1:n+1)
for(j=1:m+1)
p(i,j)=exp(x(j)+t(i));
e(i,j)=abs(u(i,j)-p(i,j));
end
end
在命令窗口中输入：
[p u e x t]=pwxywxq(0.1,0.005,10,200);
>> surf(x,t,u)
>> shading interp;
>> xlabel('x');ylabel('t');zlabel('u');
>> title('一维抛物线方程向前欧拉法数值解');
surf(x,t,p)
>> shading interp;
xlabel('x');ylabel('t');zlabel('p');
>> title('一维抛物线方程向前欧拉法精确解')
同理：
plot(x,u)
>> xlabel('x');ylabel('u');
>> title('固定时间改变x u与x 的关系数值解')
[p u e x t]=pwxywxq(0.1,0.01,10,100);
surf(x,t,u)
Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Not rendering
Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Not rendering
Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Not rendering
Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Not rendering
>> surf(x,t,e)
Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Not rendering
Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Not rendering
Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Not rendering
>>
所以空间步长与时间步长需要满足上面所说的关系
继续减小时间步长
[p u e x t]=pwxywxq(0.1,0.001,10,1000)
此为欧拉向前差分法，向后差分法请参看下一篇文章
：一维抛物线偏微分方程数值解法（2）（附matlab程序及图片）
我近期在做这个，有兴趣可以一起学习
百度账号：草随风逝。