中心差分解两点边值问题

1、一题目用中心差分格式计算如下两点边值问题已知其精确解为二理论作为模型，考虑两点边值问题：（1.1）（1.2）假定是给定的常数。1. 建立差分格式（1）.区域网格剖分首先取个节点：将区间分成个小区间：于是得到区间的一个网格剖分。记，称为网格最大步长。用表示网格内点，的集合，表示内点和界点的集合。取相邻节点的中点，称为半整数点。则由节点又构成的一个网格剖分，称为对偶剖分。（2）.微分方程的离散，建立相应差分格式用差商代替微商，将方程（1.1）在内点离散化.注意对充分光滑的，由Taylor展式有（1.3） （1.5）由（1.5）减（1.4），并除以，得（1.6）令则由（1.3）（1.6）知，边值问题

2、的解满足方程： （1.7） 其中 （1.8）为差分算子的截断误差，舍去，便得逼近边值问题（1.1）（1.2）的差分方程： （1.9）i=1,2,N-1，由方程（1.7）（1.9），截断误差可表示为 （1.10）当网格均匀，即时差分方程（1.9）简化为 （1.11）这相当于用一阶中心差商，二阶中心差商依次代替（1.1）的一阶微商和二阶微商的结果。这个方程就是中心差分格式。截断误差为： （1.12）所以截断误差按或的阶为。在本题中， ，因为r=0方程（1.11）的系数对角矩阵是三对角矩阵。我们可以用消元法或迭代法求解方程组（1.1）（1.2）式（1.11）用方程组展开：写成矩阵形式为：2.收敛性分

3、析根据（1.10）我们引进误差则误差函数满足下列差分方程：于是收敛性及收敛速度的估计问题，就归结到通过右端（截断误差）估计误差函数的问题。由（1.12）我们知，有从而差分方程满足相容条件。若引进记号，设则可将（1.9）改写为将差分解表成 （2.1）其中满足 （2.2）而满足 （2.3）先估计，由 （2.4）据差分格林公式再利用柯西不等式，有常数使 （2.5）将不等式（2.6）用于（2.5）右端，则 （2.6）解差分方程（2.2，易得）从而这样， （2.7）利用范数，从（2.7）推出 （2.8）因为因此 （2.9）联结（2.1）（2.7）及（2.9）即得差分解的先验估计： （2.10）其中不等式

4、（2.10）说明差分解连续依赖于右端和边值，因此差分格式（1.11）关于右端及边值稳定.根据定理1.1 ： 若边值问题的解u充分光滑，差分方程按满足相容条件且关于右端稳定，则差分解按收敛到边值问题的解，且有和相同的收敛阶。所以差分方程的解的收敛速度为。三程序代码：clcclfclfsyms x;a=1; %区间界点b=2; %区间界点p=exp(x); %这是p函数q=sin(x)+1+x; %这是q函数f=-exp(x)*(2*x+1)+(sin(x)+1+x)*x*(x-1);%这是f函数r=0; %这是r函数.N=10; %将区间划分的等分,这里控制！h=(b-a)/N; %这里确定步长

5、value_of_f=zeros(N-1,1);%这是fdiag_0=zeros(N-1,1);%确定A的对角元diag_1=zeros(N-2,1);%确定A的偏离对角的上对角元diag_2=zeros(N-2,1);%确定A的偏离对角的下对角元X=a:h:b;u_a=0; %边界条件u_b=2; %边界条件for j=2:N diag_0(j-1)=(subs(p,x,(X(j+1)+X(j)/2)+(subs(p,x,(X(j-1)+X(j)/2)/(h2)+(subs(q,x,X(j);end %获取对角元素for j=3:N diag_2(j-2)=-(subs(p,x,(X(j-1

6、)+X(j)/2)/(h2)-subs(r,x,X(j)/(2*h);end %获取A的第三条对角for j=2:N-1 diag_1(j-1)=-(subs(p,x,(X(j+1)+X(j)/2)/(h2)+subs(r,x,X(j)/(2*h);end %获取A的第二条对角for j=2:N; value_of_f(j-1)=subs(f,x,X(j);end %获取F值value_of_f(1)=value_of_f(1)+u_a*(subs(p,x,(X(2)+X(1)/2)/(h2);value_of_f(N-1)=value_of_f(N-1)+u_b*(subs(p,x,(X(N

7、)+X(N+1)/2)/(h2);A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1);%组装系数矩阵format longU=inv(A)*value_of_f %差分解%fprintf('%11.5f',U)fprintf('n');dx=X(2:N);precise_value=dx.*(dx-1) %精确解%fprintf('%11.5f',precise_value)deta=U-precise_value' ; %误差deta_max=max(abs(deta);%最大误差fprintf(

8、'最大的误差是%fn',deta_max)plot(X(2:N),U,'b*',X(2:N),precise_value,'r-') %差分解与精确解对比表figure();plot(X(2:N),deta) %误差图结果： X的值步长h2.12.22.32.42.52.62.72.82.9最大误差0.10.110150.240260.390330.560370.750380.960381.190311.440221.710120.0003800.050.110030.240060.390080.560090.750090.960081.190071.440051.710030.0000950.0250.110000.240010.390020.560020.750020.960021.190021.440011.710010.0000240.01250.110000.240000.390000.560010.750010.960001.190001.440001.710000.000006精确解0.110000.240000.390000.560000.750000.960001.190001.440001.710000以下仅给出步长为N=20,h=0.05的精确值和差分值图与误差图，其