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

一维抛物线偏微分方程数值解法（1） 解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法） Ut-Uxx=0, 0x1,0t=1（Ut-aUxx=f(x,t),a0) U(x,0)=e^x, 0=x=1, U(0,t)=e^t,U(1,t)=e^(1+t), 0t=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),a0) %不用解线性方程组，由下一层(时间层）的值就直接得到上一层的值 %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),a0) %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(kkmax) 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