2022年一维抛物线偏微分方程数值解法.docx
名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -一维抛物线偏微分方程数值解法（3）上一篇参看 一维抛物线偏微分方程数值解法（2）（附图及 matlab程序）解一维抛物线型方程 理论书籍可以参看孙志忠：偏微分方程数值解法）Ut-Uxx=0, 0<x<1,0<t<=1（Ut-aUxx=fx,t,a>0 Ux,0=ex, 0<=x<=1, U0,t=et,U1,t=e1+t, 0<t<=1 精确解为： Ux,t=ex+t; 此种方法精度为 oh12+h22 一：用追逐法解线性方程组（仍可以用迭代法解）Matlab 程序 function u p e x t=CNh1,h2,m,n %Crank-Nicolson 格式差分法解一维抛物线型偏微分方程%此程序用的是追逐法解线性方程组 %h1为空间步长, h2 为时间步长 %m,n 分别为空间,时间网格数 %p为精确解, u 为数值解 ,e 为误差x=0:m*h1+0; x0=0:m*h1; %定义 x0,t0是为了 fx,t=0的情形 % t=0:n*h2+0; t0=0:n*h2+1/2*h2; syms f ;fori=1:n+1%fi,j=fx0j,t0i=0%forj=1:m+1 fi,j=0; end end for i=1:n+1 ui,1=expti; ui,m+1=exp1+ti; end for i=1:m+1 u1,i=expxi; end r=h2/h1*h1;fori=1:n %外循环,先固定每一时间层,每一时间层上解一线性方程组%.a1=0;b1=1+r;c1=-r/2;d1=r/2*ui+1,1+ui,1+h2*fi,j +1-r*ui,2+r/2*ui,3; for k=2:m-2ak=-r/2;bk=1+r;ck=-r/2;dk=h2*fi,j+r/2*ui,k+1-r . *ui,k+1+r/2*ui,k+2;细心整理归纳 精选学习资料 %输入部分系数矩阵,为0的矩阵元素不输入% 第 1 页,共 9 页 - - - - - - - - - - - - - - - - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -endam-1=-r/2;bm-1=1+r;dm-1=h2*fi,j+r/2*ui,m+1+ui+1,m+1 . +r/2*ui,m-1+1-r*ui,m;fork=1:m-2 %开头解线性方程组消元过程 ak+1=-ak+1/bk; bk+1=bk+1+ak+1*ck; dk+1=dk+1+ak+1*dk; end ui+1,m=dm-1/bm-1; %回代过程 % for k=m-2:-1:1 ui+1,k+1=dk-ck*ui+1,k+2/bk;end endfori=1:n+1%p为精确解forj=1:m+1 pi,j=expxj+ti; ei,j=absui,j-pi,j;%e为误差end endu p e x t=CN0.1,0.005,10,200;surfx,t,e; shading interp; >> xlabel'x'ylabel't'zlabel'e' >> title' 误差曲面 ' 细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 2 页,共 9 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -plotx,e plott,e 细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 3 页,共 9 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -误差较向前欧拉法减小一半但是运行时间较长,约39 秒,而前两次运行只需l 秒左右；,运算量也提高五倍 第 4 页,共 9 页 u p e x t=CN0.01,0.01,100,100; 运行需三分钟左右,误差比前次提高五倍细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -u p e x t=CN0.1,0.1,10,10;surfx,t,e u p e x t=CN0.1,0.2,10,5;surfx,t,e 运行需要 2秒；精度仍是挺高的；细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 5 页,共 9 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -误差仍可以接受 此种方法精度高,运算量较大二：用迭代法解线性方程组：Matlab 程序如下：function u e p x t k=CN1h1,h2,m,n,kmax,ep % 解抛物线型一维方程 C-N 格式 （ Ut-aUxx=fx,t,a>0 %用g-s 高斯 - 赛德尔 迭代法解 %kmax为最大迭代次数%m,n 为x,t方向的网格数,例如（2-0 ）/0.01=200;%e为误差, p为精确解 syms temp ; u=zerosn+1,m+1; x=0+0:m*h1; t=0+0:n*h2; for i=1:n+1 ui,1=expti; ui,m+1=exp1+ti; end for i=1:m+1 u1,i=expxi; endfori=1:n+1 第 6 页,共 9 页 细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -for j=1:m+1 fi,j=0; end end a=zerosn,m-1;r=h2/h1*h1; %此处 r=a*h2/h1*h1）；a=1.fork=1:kmaxfori=1:nforj=2:m temp=r/2*ui,j-1+1-r*ui,j+r/2*ui,j+1+h2*fi,j+r/2*ui+1,j-1+r/2*ui+1,j+1/1+r; ai+1,j=temp-ui+1,j*temp-ui+1,j; ui+1,j=temp;%此处留意是 ui+1,j,而不是 ui+1,j+1%end end ai+1,j=sqrtai+1,j;ifk>kmax;breakend if maxmaxa<ep break ; end end for i=1:n+1 for j=1:m+1 pi,j=expxj+ti; ei,j=absui,j-pi,j; end end u e p x t k=CN10.1,0.005,10,200,10000,1e-10; 运行速度： 1秒 迭代次数 k = 81 surfx,t,e 细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 7 页,共 9 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -其次幅图为三角追逐法解方程作出的图,两者几乎一样；细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 8 页,共 9 页 - - - - - - - - - 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -由于迭代法速度很快,所以可以将区间分得更小u e p x t k=CN10.01,0.01,100,100,10000,1e-12;surfx,t,e;shading interp; k=6903细心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 9 页,共 9 页 - - - - - - - - -
