偏微分中心差分格式实验报告(含matlab程序)

二阶常微分方程的中心差分求解学校:中国石油大学（华东）理学院姓名：张道德一、实验目的1、构造二阶常微分边值问题:22,(),(),duLuqufaxbdxuaub其中,qf为[,]ab上的连续函数，0;,q为给定常数的中心差分格式；2、根据中心差分格式求解出特定例题的数值解，并与该例题的精确解进行比较。二、中心差分格式的构造将区间[,]ab分成N等分，分点为：0,1,2,,ixaihiN()/hbaN。于是我们得到区间的一个网络剖分。称为网格的节点称为步长。得中心差分格式为：11202,1,2,,1,,.iiihiiiiNuuuLuqufiNhuu其中式中(),()iiiiqqxffx。三、差分格式求解根据中心差分格式可以构造出：1112222222233322212211210012101201001200NNNufqhhufqhhhufqhhhqufhh可以看出系数矩阵为三对角矩阵，而对于系数矩阵为三对角矩阵的方程组可以用“追赶法”求解，则可以得出二阶常微分方程问题的数值解。四、举例求解我们选取的二阶常微分方程边值问题为：222242,01(0)1,(1),xduLuxuexdxuue其精确解为：2xue。则我们具体求解出的解如下：1、当10N时，数值解与精确解为：(1)表1、10N时，数值解与精确解统计表x的值0.10.20.30.40.5u的数值解1.0110691.0427441.0969041.1768961.28789u的精确解1.010051.0408111.0941741.1735111.284025两者之差0.0010190.0019340.0027290.0033850.003864x的值0.60.70.80.9u的数值解1.4374431.6363631.9000012.250209u的精确解1.4333291.6323161.8964812.247908两者之差0.0041140.0040460.003520.002301将两者绘于同一图中如下：(2)结论：可以看出数值解与精确解之间的误差很小，在210这样一个数量级上。我们也可以求出的|-|数值解精确解范数来，如下：Norm1(|-|数值解精确解)=0.0269;Norm2(|-|数值解精确解)=0.0095;Normoo(|-|数值解精确解)=0.0041;所以，可以得出中心差分格式求解该方程效果挺好。2、当20N时，将数值解与精确解绘于同一图像中，如下：0.10.20.30.40.50.60.70.80.911.522.5x的值u的数值解与精确解图1、N=10时，数值解与精确解图像数值解精确解(2)结论：可以看出数值解与精确解之间的误差很小，在410这样一个数量级上。我们也可以求出的|-|数值解精确解范数来，如下：Norm1(|-|数值解精确解)=0.0027;Norm2(|-|数值解精确解)=43.018810;Normoo(|-|数值解精确解)=54.166910;所以，可以得出中心差分格式求解该方程效果挺好。五、程序程序100.10.20.30.40.50.60.70.80.9111.21.41.61.822.22.42.62.8x的值u的数值解与精确解图2、当N=100时，u的数值解与精确解数值解精确解%**************************************************************%f221.mfunction[q,f]=f211(x)%q函数，f函数q=4*x^2;f=-2*exp(x^2);%***************************************************************程序2%********************************************************************%追赶法function[x]=zhuiganfa(a,b,c,d)%对角线下方的元素，个数比A少一个%%对角线元素%对角线上方的元素，个数比A少一个%d为方程常数项%用追赶法解三对角矩阵方程r=size(a);m=r(2);r=size(b);n=r(2);ifsize(a)~=size(c)|m~=n-1|size(b)~=size(d)error('变量不匹配,检查变量输入情况!');end%%%LU分解u(1)=b(1);fori=2:nl(i-1)=a(i-1)/u(i-1);u(i)=b(i)-l(i-1)*c(i-1);v(i-1)=(b(i)-u(i))/l(i-1);end%追赶法实现%%%求解Ly=d，追的过程y(1)=d(1);fori=2:ny(i)=d(i)-l(i-1)*y(i-1);end%%%求解Ux=y,赶的过程x(n)=y(n)/u(n);fori=n-1:-1:1x(i)=y(i)/u(i);x(i)=(y(i)-c(i)*x(i+1))/u(i);end%********************************************************************程序3%*********************************************************************%ODE2.mfunction[x]=ODE2(x0,xN,u0,uN,N)%中心差分求解%x0,xN初始条件%u0,uN边值条件%N等分数%步长h=1/N;%%a(1:N-2)=-1/h^2;%对角线下方的元素，个数比A少一个fori=1:N-1z(i)=x0+i*h;[q,f]=f211(z(i));b(i)=2/h^2+q;%对角线元素d(i)=f;end%对对角线元素进行调整d(1)=d(1)+u0/h^2;d(N-1)=d(N-1)+uN/h^2;%%[x]=zhuiganfa(a,b,a,d);%数值解%**********************************************************************程序4%**************************************************************%main_chapter.m%x0,xN初始条件%u0,uN边值条件%N等分数%h步长x0=0;xN=1;u0=1;uN=exp(1);N=10;h=1/N;[x]=ODE2(x0,xN,u0,uN,N);%数值解z=x0+h:h:xN-h;y=exp(z.^2);%真解nr1(x-y)plot(z,x,'og',z,y,'.r')%***************************************************************