常微分方程模型及其数值解2008.ppt

常微分方程数值解0、导言另一个是方程的求解和结果分析。对一些常系数的或特殊函数形式的微分方程，往往能得到解析解，这对实际问题的分析和应用都是有利的,但是大多数变系数的、非线性函数形式的微分方程都是求不出解析解的,此时就需要应用求解微分方程的另一个重要方法──数值解法。
本章简要介绍有关微分方程模型的概念，微分方程的数值解法和图解法，主要介绍若干建模实例，通过它们展示微分方程模型的建模步骤及解决实际问题的全过程。欧拉方法和龙格—库塔方法所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值.对右边的积分应用左矩形公式，则有对右边的积分应用左矩形公式，则有利用Euler方法求初值问题分别取步长h=0.2,0.1,0.05,计算结果如下function[outx,outy]=MyEuler(fun,x0,xt,y0,PointNum)
%MyEuler用前向差分的欧拉方法解微分方程
%fun表示f（x，y）
%x0,xt表示自变量的初值和终值
%y0表示函数在x0处的值,其可以为向量形式
%PointNum表示自变量在[x0,xt]上取的点数
ifnargin<5|PointNum<=0%如果函数仅输入4个参数值，则PointNum默认值为100
PointNum=100;
end
ifnargin<4%y0默认值为0
y0=0;
endh=(xt-x0)/PointNum;%计算步长h
x=x0+[0:PointNum]'*h;%自变量数组
y(1,:)=y0(:)';%将输入存为行向量，输入为列向量形式
fork=1:PointNum
f=feval(fun,x(k),y(k,:));%计算f(x,y)在每个迭代点的值
f=f(:)';
y(k+1,:)=y(k,:)+h*f;%对于所取的点x迭代计算y值
end
outy=y;
outx=x;
%plot(x,y)%画出方程解的函数图
functionf=myfun011(x,y)
f=1/(1+x*x)-2*y*y;hEuler中点公式则不然,计算yn+1时需用到前两步的值yn,yn-1,称其为两步方法,两步以上的方法统称为多步法.从数值积分的角度来看,梯形公式由迭代法收敛的角度看，当y=y-2x/y,0x1计算结果如下:function[Xout,Yout]=MyEulerPro(fun,x0,xt,y0,PointNumber)
%MyEulerPro用改进的欧拉法解微分方程
ifnargin<5|PointNumber<=0%如果函数仅输入4个参数值，则PointNumer默认值为100
PointNumer=100;
end
ifnargin<4%y0默认值为0
y0=0;
endh=(xt-x0)/PointNumber;%计算所取的两离散点之间的距离
x=x0+[0:PointNumber]'*h;%表示出离散的自变量x
y(1,:)=y0(:)';
fori=1:PointNumber%迭代计算过程
f1=h*feval(fun,x(i),y(i,:));
f1=f1(:)';
f2=h*feval(fun,x(i+1),y(i,:)+f1);
f2=f2(:)';
y(i+1,:)=y(i,:)+1/2*(f1+f2);
end
Xout=x;
Yout=y;functionf=myfun012(x,y)
f=y-2*x/y;n在节点xn+1的误差y(xn+1)-yn+1,不仅与yn+1这一步计算有关,而且与前n步计算值yn,yn-1,…,y1都有关.另外,在yn=y(xn)的条件下,考虑到y(x)=(x,y(x)),则有yn+1=yn+h(xn,yn)所以,改进的Euler方法是二阶方法.称之为p阶Taylor展开方法.所以,此差分公式是p阶方法.构造差分公式由于一般地,参数由(8)确定的一族差分公式(7)统称为二阶R-K方法.四阶标准R-K公式解四阶标准R-K公式为function[x,y]=MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)
%Runge-Kutta方法解微分方程形为y’(t)=f(x,y(x))
%此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式
%x范围为[x0,xt],初值为y0,PointNum为离散点数,varargin为可选输入项可传适当参数给函数f(x,y)
ifnargin<4|PointNum<=0
PointNum=100;
end
ifnargin<3
y0=0;
endy(1,:)=y0(:)';%初值存为行向量形式
h=(xt-x0)/(PointNum-1);%计算步长
x=x0+[0:PointNum]'*h;%