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

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

二阶常微分方程的中心差分求解学校:中国石油大学(华东)理学院 姓名：张道德实验目的构造二阶常微分边值问题:其中为上的连续函数，为给定常数的中心差分格式；根据中心差分格式求解出特定例题的数值解，并与该例题的精确解进行比较。中心差分格式的构造将区间分成等分，分点为：。于是我们得到区间的一个网络剖分。称为网格的节点称为步长。得中心差分格式为：其中式中。差分格式求解根据中心差分格式可以构造出：可以看出系数矩阵为三对角矩阵，而对于系数矩阵为三对角矩阵的方程组可以用“追赶法”求解，则可以得出二阶常微分方程问题的数值解。举例求解我们选取的二阶常微分方程边值问题为： 其精确解为：。则我们具体求解出的解如下：当时，数值解与精确解为：(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将两者绘于同一图中如下： (2)结论：可以看出数值解与精确解之间的误差很小，在这样一个数量级上。我们也可以求出的范数来，如下：Norm1()=0.0269;Norm2()=0.0095;Normoo()=0.0041;所以，可以得出中心差分格式求解该方程效果挺好。2、当时，将数值解与精确解绘于同一图像中，如下：(2)结论：可以看出数值解与精确解之间的误差很小，在这样一个数量级上。我们也可以求出的范数来，如下：Norm1()=0.0027;Norm2()=;Normoo()=;所以，可以得出中心差分格式求解该方程效果挺好。程序程序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);if size(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]=zhui