微分方程的matlab求解.

主页下一页上一页数学实验定义：含有导数的方程称为微分方程。如f(x,y(x),y’(x))=0微分方程模型1、微分方程的一般形式：F(x,y,y’,…,y(n))=0隐式或y(n)=f(x,y,y’,…,y(n-1))显式主页下一页上一页数学实验特殊情形：d(,)dyftyt2、一阶微分方程组的一般形式：));(,),(),(,());(,),(),(,(212111xyxyxyxfdxdyxyxyxyxfdxdymmmm初始条件：y(x0)=y0微分方程模型主页下一页上一页数学实验③图形解tyo①简单的微分方程。②复杂、大型的微分方程。返回①解析解y=f(t)②数值解(ti,yi)•欧拉方法•改进欧拉方法•梯形法•龙格-库塔法微分方程求解方法简介主页下一页上一页数学实验微分方程数值解1、欧拉法2、龙格—库塔法数值求解思想：(变量离散化)引入自变量点列{xn}→{yn}，在x0x1x2…xn…上求y(xn)的近似值yn.通常取等步长h，即xn=x0+n×h，或xn=xn-1+h，（n=1,2,…)。00)(),,(yxyyxfdxdy主页下一页上一页数学实验1)向前欧拉公式:（y’=f(x,y)）y(xn+1)y(xn)+hf(xn,y(xn))(迭代式)yn+1yn+hf(xn,yn)(近似式)特点：f(x,y)取值于区间[xn,xn+1]的左端点.')()(1yhxyxynn在小区间[xn,xn+1]上用差商代替微商(近似),1、欧拉方法主页下一页上一页数学实验yn+1yn+hf(xn+1,yn+1)特点：①f(x,y)取值于区间[xn,xn+1]的右端点.②非线性方程,称‘隐式公式’。yn+1=yn+hf(xn,yn)2)向后欧拉公式方法：迭代（y’=f(x,y)）x=[];y=[];x(1)=x0;y(1)=y0;forn=1:kx(n+1)=x(n)+n*h;y(n+1)=y(n)+h*f(x(n),y(n));(向前）end1、欧拉方法主页下一页上一页数学实验1),(1.0,1)0(,1'xyyxfhyxyy其中例1观察向前欧拉、向后欧拉算法计算情况。与精确解进行比较。误差有多大？解：1）解析解：y=x+e-x1、欧拉方法主页下一页上一页数学实验2）向前欧拉法：yn+1=yn+h(-yn+xn+1)=(1-h)yn+hxn+h3）向后欧拉法：yn+1=yn+h(-yn+1+xn+1+1)转化yn+1=(yn+hxn+1+h)/(1+h)y’=f(x,y)=-y+x+1;1、欧拉方法主页下一页上一页数学实验x1(1)=0;y1(1)=1;y2(1)=1;h=0.1;(died.m)fork=1:10x1(k+1)=x1(k)+h;y1(k+1)=(1-h)*y1(k)+h*x1(k)+h;y2(k+1)=(y2(k)+h*x1(k+1)+h)/(1+h);endx1,y1,y2,（y1——向前欧拉解，y2——向后欧拉解）x=0:0.1:1;y=x+exp(-x)（解析解）plot(x,y,x1,y1,'k:',x1,y2,'r--')1、欧拉方法主页下一页上一页数学实验x精确解向前欧拉向后欧拉01110.11.004811.00910.21.01871.011.02640.31.04081.0291.05130.41.07031.05611.08300.51.10651.09051.12090.61.14881.13141.16450.71.19661.17831.21320.81.24931.23051.26650.91.30661.28741.324111.36791.34871.3855（1）步长h=0.1的数值解比较表计算结果主页下一页上一页数学实验（2）步长h=0.01的数值解比较表x精确解向前欧拉向后欧拉01110.11.00481.00441.00530.21.01871.01791.01950.31.04081.03971.04190.41.07031.06901.07170.51.10651.10501.10800.61.14881.14721.15040.71.19661.19481.19830.81.24931.24751.25110.91.30661.30471.308411.36791.36601.3697结论：显然迭代步长h的选取对精度有影响。主页下一页上一页数学实验00.5111.11.21.31.4xx+exp(-x)图形显示Z有什么方法可以使精度提高？向后欧拉法返回主页下一页上一页数学实验梯形公式,...2,1,0)],(),([2111nyxfyxfhyynnnnnn改进欧拉公式),(),()(21121211hkyxfkyxfkkkhyynnnnnnyn+hf(xn,yn)返回主页下一页上一页数学实验x精确解向前欧拉向后欧拉改进欧拉011110.11.004811.00911.0050.21.01871.011.02641.0190.31.04081.0291.05131.04120.41.07031.05611.08301.07080.51.10651.09051.12091.10710.61.14881.13141.16451.14940.71.19661.17831.21321.19720.81.24931.23051.26651.25000.91.30661.28741.32411.307211.36791.34871.38551.3685步长h=0.1的数值解比较表使用改进欧拉公式的例主页下一页上一页数学实验2、龙格-库塔法龙格-库塔法是利用泰勒展式将y(x+h)在x处展开，并取其前面若干项来近似y(x+h)而得到公式y(x+h)y(x)+hj(x,y(x),h)如果y(xn)yn，则y(xn+1)的近似值为：yn+1=yn+hj(xn,yn,h),n=0,1,…若y(x+h)[y(x)+hj(x,y(x),h)]=O(hp+1)，则称以上迭代公式为p阶公式，p的大小反映了截断误差的高低，高阶高精度。要得到一个p阶公式，关键在于如何选取j(x,y(x),h)使之满足阶的要求。返回主页下一页上一页数学实验微分方程图解法欲将微分方程解的全局信息形象化、直观化。对于一阶微分方程dy/dx=f(x,y)，如果给出平面上任意一点(x,y)，就能够确定出解y=f(x)在该点(x,y)处的斜率f(x,y)。从图象上看，给出平面上的一系列点，通过每一点(x0,y0)，可以画出一条通过点(x0,y0)、斜率为f(x0,y0)的短直线。这样的短直线布满整个坐标平面，形成的图形就称为斜率场或方向场。主页下一页上一页数学实验-10-8-6-4-20246810-10-8-6-4-20246810xy斜率场:dy/dx=sinx*siny微分方程图解法主页下一页上一页数学实验相平面轨迹表示微分方程的解22220000()()()()(),()dxtyxxxydtdytyxyxydtxtxyty微分方程图解法利用微分方程的数值解法，可以得到其数值解：(x(t)，y(t))在t取离散值时的取值列向量X,Y；然后分别独立地作出函数x(t)和y(t)的曲线，如图4.2，其初值条件为(5,5)。主页下一页上一页数学实验微分方程图解法00.511.522.533.544.55-20246tx数值解曲线00.511.522.533.544.55-20246ty主页下一页上一页数学实验微分方程图解法如果撇开自变量的取值T，直接利用X,Y的分量作为坐标，就可以在xoy平面上画出解的轨迹，称为相平面轨迹图。-2-1012345-2-1012345xy相平面图主页下一页上一页数学实验微分方程图解法-2-1012345-2-1012345斜率场xy主页下一页上一页数学实验MATLAB软件求解解析解y=dsolve('eqn1','eqn2',…,'c1',…,'x')微分方程组初值条件指明变量注意：①y'Dy，y''D2y②自变量名可以省略，默认变量名‘t’。主页下一页上一页数学实验例①1)0(,12yydxdy输入：y=dsolve('Dy=1+y^2')y1=dsolve('Dy=1+y^2','y(0)=1','x')输出：y=tan(t-C1)（通解）y1=tan(x+1/4*pi)（特解）MATLAB软件求解主页下一页上一页数学实验例②常系数的二阶微分方程0)0(',1)0(,03'2''yyyyyy=dsolve('D2y-2*Dy-3*y=0','x')y=dsolve('D2y-2*Dy-3*y=0','y(0)=1,Dy(0)=0','x')输入：y=C1*exp(-x)+C2*exp(3*x)y=3/4*exp(-x)+1/4*exp(3*x)结果：主页下一页上一页数学实验x=dsolve('D2x-(1-x^2)*Dx+x=0','x(0)=3,Dx(0)=0')例③非常系数的二阶微分方程0)0(',3)0(,0)()('))(1()(''2xxtxtxtxtx无解析表达式！主页下一页上一页数学实验x=dsolve('(Dx)^2+x^2=1','x(0)=0')例④非线性微分方程0)0(,1)()('22xtxtxx=sin(t)-sin(t)若欲求解的某个数值解，如何求解？t=pi/2;eval(x)MATLAB软件求解主页下一页上一页数学实验输入：[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y')[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y','x(0)=0,y(0)=1')例④1)0(0)0(3443yxyxdtdyyxdtdx输出：x=1/2*exp(7*t)-1/2*exp(-t)y=1/2*exp(-t)+1/2*exp(7*t)返回MATLAB软件求解主页下一页上一页数学实验MATLAB软件求解数值解其中（1）Fun表示由微分方程(组)写成的m文件名；（2）y0表示为函数的初值；（3）options用于设定误差限(可以默认)。程序为options=odeset('reltol',rt,'abstol',at)这里的rt和at分别为设定的相对误差和绝对误差。(rt=1e-7)[t,y]=ode45('Fun',[t0，tf],y0,options)主页下一页上一页数学实验1）首先建立M-文件（weif.m）functionf=weif(x,y)f=-y+x+1;2）求解：[x，y]=ode23(‘weif’,[0,1],1)3)作图形:plot(x,y,‘r’);4)与精确解进行比较holdonezplot(‘x+exp(-x)’,[0,1])例1y’=-y+x+1，y(0)=1标准形式：y’=f(x,y)范例主页下一页上一页数学实验使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意:()(1)(1)(,,,,)(0),(0),,(0)nnnyftyyyyyy选择一组状态变量(1)12,,,nnxyxyxy122312,,(,,,,)nnxxxxxftxxx主页下一页上一页数学实验注意1、建立M文件函数functionxdot=fun(t,x,y)xdot=[x2(t)；x3(t)；…；f(t,x1(t),x2(t),…xn(t))]；2、数值计算（执行以下命令）[t,x1,x2,…,xn]=ode45(‘fun',[t0,tf],[x1(0),x2(0),…,xn(0)])122312,,(,,,,)nnxxxxxftxxx主页下一页上一页数学实验例2Vanderpol方程：0)0(',3)0(0)()('))(1()(''2