试射法&差分法

----试射法和差分法

摘要

本文运用边值问题的两种解法，试射法和差分法，解二阶微分方程。以一个实例进行matlab编程实现，对这两种解法进行进一步的了解和比较。

关键词：常微分方程边值问题，试射法，差分法

一、问题分析

   对于二阶微分方程

边界条件可分为三类：①            

                    ②            

                    ③    

                      （，）

分别称为第一、二、三类边界条件。不同边界条件对用的微分方程求解问题分别称为第一、二、三边值问题。本文讨论第一边值问题。边值问题解的存在性不仅与方程有关，也与边界条件有关，所以我们始终假定边值问题有解，且解具有所需的光滑性。

1.1试射法

1、基本思想：

    将第一边值问题

的求解转化为初值问题的求解。但初始条件无法确定，因此考虑用一组参数逐步逼近，即考虑初值问题使其解满足。这样，对给定的允许误差，当或时，即可作为所求之近似解。

2、实用方法选择：

    （1）根据对边值问题性质的分析及问题的实际背景，对提供一个初始估计值，例如可取

然后解初值问题，得解，记。

（2）取使满足即取，再解初值问题，得，记。若，则为所求解，否则进行下一步。

（3）利用线性插值法由、求得，再取进行试射，得新解，记。若，则为所求解；否则令，，，重复此步计算过程，直至达到目的为止。

1.2差分法

1、基本思想：

接编制问题的差分方法是一种常用方法，用差商代导数，将微分方程离散化为差分方程，然后用适当的方法求解。

2、具体步骤：

将求解区间分成等分，步长，分点称为节点，在内节点上将导数，用差商表示

在点列出方程

代入上式整理得到近似方程

其中，，。

此外，由边界条件得

将两式联立在一起恰好构成个方程、个未知量的线性方程组，称为边值问题的差分方程，其解称为边值问题的差分解，称为阶段误差。整理两式可得三对角方程组

其中

此方程组可用追赶法求解。

二、算法实现及结果

例题  用有限差分法解例题2.1中的边值问题：

2.1试射法

利用RK方法计算二阶线性微分方程的边值问题，由微分方程知识可以知道，对于二阶微分方程可以化为一阶微分方程组，func1为第一个微分方程组第二个函数，func2为第二个微分方程组第二个函数，初始时刻为a，结束时为b，初始端点函数值为ya，末端点函数值yb，N为分成的区间数目。

用matlab实现算法的程序：

function [x,y]=lsh(func1,func2,a,b,ya,yb,N)

h=(b-a)/N;

u(1,1)=ya;

u(1,2)=0;

v(1,1)=0;

v(1,2)=1;

for i=1:N

    x(i,:)=a+(i-1)*h;

    K1=h*u(i,2);

    L1=h*feval(func1,x(i,:),u(i,1),u(i,2));

    K2=h*(u(i,2)+1/2*L1);

    L2=h*feval(func1,x(i,:)+1/2*h,u(i,1)+1/2*K1,u(i,2)+1/2*L1);

    K3=h*(u(i,2)+1/2*L2);

    L3=h*feval(func1,x(i,:)+1/2*h,u(i,1)+1/2*K2,u(i,2)+1/2*L2);

    K4=h*(u(i,2)+L3);

    L4=h*feval(func1,x(i,:)+h,u(i,1)+K3,u(i,2)+L3);

    u(i+1,1)=u(i,1)+1/6*(K1+2*K2+2*K3+K4);

    u(i+1,2)=u(i,2)+1/6*(L1+2*L2+2*L3+L4);

    k1=h*v(i,2);

    l1=h*feval(func2,x(i,:),v(i,1),v(i,2));

    k2=h*(v(i,2)+1/2*l1);

    l2=h*feval(func2,x(i,:)+1/2*h,v(i,1)+1/2*k1,v(i,2)+1/2*l1);

    k3=h*(v(i,2)+1/2*l2);

    l3=h*feval(func2,x(i,:)+1/2*h,v(i,1)+1/2*k2,v(i,2)+1/2*l2);

    k4=h*(v(i,2)+l3);

    l4=h*feval(func2,x(i,:)+h,v(i,1)+k3,v(i,2)+l3);

    v(i+1,1)=v(i,1)+1/6*(k1+2*k2+2*k3+k4);

    v(i+1,2)=v(i,2)+1/6*(l1+2*l2+2*l3+l4);    

end

u

v

x(N+1,:)=x(N,:)+h;

y(1,1)=ya;

y(1,2)=(yb-u(N+1,1))/v(N+1,1);

当，时，计算结果如下表所示。表中列出的值逼近的值，逼近的值，且值逼近

，

表示逐点绝对误差。

二阶线性微分方程边值问题有限差分算法，y''=py'+qy+r，y(a)=alpha,y(b)=beta，输入参数，p,q,r 微分方程描述参数， a,b 微分区间[a,b]，alpha,beta 上下限对应的边值，n 等分区间个数，输出参数，t 区间点坐标，y 函数值。

用matlab实现算法的程序：

function [t,y]=findiff(p,q,r,a,b,alpha,beta,n)

vb=zeros(1,n-1);%等式右边

va=zeros(1,n-2);% +1对角线

vc=zeros(1,n-2);% -1对角线

vd=zeros(1,n-1);% 主对角线

h=(b-a)/n;

vt=a+h:h:a+h*(n-1);

vb=-h^2*r(vt);

vb(1)=vb(1)+(1+h/2*p(vt(1)))*alpha;

vb(n-1)=vb(n-1)+(1-h/2*p(vt(n-1)))*beta;% 计算主对角线

vd=2+h^2*q(vt);% 计算+1对角线

vta=vt(1,2:n-1);

va=-1-h/2*p(vta);% 计算-1对角线

vtc=vt(1,1:n-2);

vc=-1+h/2*p(vtc);

x=trisolve(va,vd,vc,vb);

t=[a,vt,b]';

y=[alpha x beta]';

function x=trisolve(va,vd,vc,vb)

%追赶法求解线性三对角方程组Ax=b

%va +1对角线

%vd 主对角线

%vc -1对角线

%vb 方程右边

n=length(vb);

for kk=2:n

    mult=va(kk-1)/vd(kk-1);

    vd(kk)=vd(kk)-mult*vc(kk-1);

    vb(kk)=vb(kk)-mult*vb(kk-1);

end

x(n)=vb(n)/vd(n);

for kk=n-1:-1:1

    x(kk)=(vb(kk)-vc(kk)*x(kk+1))/vd(kk);

end

取，。值逼近，表示逐点绝对误差。

本文运用试射法和差分法解二阶微分方程。从例题求解的结果（中上文两表的结果及其精度）可知，在相同的步长下试射法的精度高于差分法，其原因是：试射法将边值问题的求解转化为响应的初值问题，使用了具有较高精度的Runge-Kutta方法。由于Runge-Kutta算法具有的精度，而有限差分法的局部截断误差阶仅有，此情况下试射法的精度高于差分法。

四、参考文献

[1]曹德欣 曹璎珞，《计算方法》（第二版），中国矿业大学出版社，2001。

[2]王正盛.《Matlab与科学计算》.国防工业出版社，2011.8。