1、 一維拋物線偏微分方程數值解法（1）解一維拋物線型方程(理論書籍可以參看孫志忠：偏微分方程數值解法）Ut-Uxx=0, 0<x<1,0<t<=1（Ut-aUxx=f(x,t),a>0) U(x,0)=ex, 0<=x<=1,U(0,t)=et,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)% 解拋物線型一維

2、方程 向前歐拉格式（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);endfor(i=1:m+1) u(1,i)=exp(x(i);endfor(i=1:n+1)for(j=1:m+1) f(i,j)=0;endendr=h2/(h1*h1);%此處r=a*h2/(h1

3、*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);endendfor(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);endend或者：function u e p x t k=paowuxianyiweixq(h1,h2,m,n,kmax,ep)% 解拋物線型一維方程 向前歐拉格式 （Ut-aUxx=f(x,t),a>0)%kmax為最大迭代次數

4、%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);endfor(i=1:m+1) u(1,i)=exp(x(i);endfor(i=1:n+1)for(j=1:m+1) f(i,j)=0;endenda=zeros(n,m-1);r=h2/(h1*h1);%此處r=a*h2/(h1*h1）；a=1 要求r<=1/2差分格式才穩定for(k=1:kma

5、x)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;endend a(i,j)=sqrt(a(i,j);if(k>kmax)break;endif(max(max(a)<ep)break; endendfor(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);endend在命令窗口中輸入： p

6、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');>> titl

7、e('一維拋物線方程 向前歐拉法 精確解')同理： 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. Notrendering Warning: Axis limits ou

8、tside float precision, use ZBuffer or Painters instead. Notrendering Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Notrendering Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Notrendering >> surf(x,t,e)Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Notrendering Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Notrendering Warning: Axis limits outside float precision, use ZBuffer or Painters instead. Notrendering >> 所以空間步長與時