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

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

二阶常微分方程的中心差分求解学校:中国石油大学(华东)理学院姓名：张道德一、实验目的1、构造二阶常微分边值问题:2,(),(),duLqfaxbxuab其中为上的连续函数，为给定常,qf[,]0;,q数的中心差分格式；2、根据中心差分格式求解出特定例题的数值解，并与该例题的精确解进行比较。二、中心差分格式的构造将区间分成等分，分点为：[,]abN0,12,ixih。于是我们得到区间的一个网络剖分。称为()/h网格的节点称为步长。得中心差分格式为：1120,12,,,.iiihiiiNuLqufNh其中式中。(),()iiiiqxfx三、差分格式求解根据中心差分格式可以构造出：11122222233322121001010NNufqhhfufqhhquf可以看出系数矩阵为三对角矩阵，而对于系数矩阵为三对角矩阵的方程组可以用“追赶法”求解，则可以得出二阶常微分方程问题的数值解。四、举例求解我们选取的二阶常微分方程边值问题为：2224,01(0)1,(),xduLxeue其精确解为：。则我们具体求解出的解如下：2x1、当时，数值解与精确解为：10N(1)表1、时，数值解与精确解统计表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将两者绘于同一图中如下：0.10.20.30.40.50.60.70.80.911.522.5x、、u、、、、、、、、、1、N=10、、、、、、、、、、、、、、、、、(2)结论：可以看出数值解与精确解之间的误差很小，在这样一个数量级上。我们也可以求出的210范数来，如下：|-|数值解精确解Norm1()=0.0269;|-|数值解精确解Norm2()=0.0095;||数值解精确解Normoo()=0.0041;|-|数值解精确解所以，可以得出中心差分格式求解该方程效果挺好。2、当时，将数值解与精确解绘于同一图像中，如20N下：00.10.20.30.40.50.60.70.80.9111.21.41.61.822.22.42.62.8x、、u、、、、、、、、、2、、N=100、、u、、、、、、、、、、、、、、(2)结论：可以看出数值解与精确解之间的误差很小，在这样一个数量级上。我们也可以求出的410范数来，如下：|-|数值解精确解Norm1()=0.0027;|-|数值解精确解Norm2()=;||数值解精确解43.018Normoo()=;|-|数值解精确解569所以，可以得出中心差分格式求解该方程效果挺好。五、程序程序1%**************************************************************%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 )%***************************************************************