matlab在常微分方程数值解中的应用

MATLAB在常微分方程数值解中的应用【摘要】许多现实问题都可以通过微分方程的形式进行表示，传统解微分方程的方法有近似分析解法、表解法和图解法，这些方法需对其进行大量的假设，而使得数学模型有一定的失真，有一定的局限性。数值解法利用计算机，使得求解更精确、效率更高，而MATLAB是一种数学软件包，有高级编程格式，使得计算结果更具有可信性，因此微分方程的求解及MATLAB在其中的应用具有实际意义。本文对常微分方程数值解问题作进一步探讨，并应用MATLAB对其中难解的改进Euler法和Runge-Kutta法进行编程实现，程序简洁、直观，求解速度快、方法实用性较强。【关键词】常微分方程数值解MATLABEuler法龙格-库塔方法ode45ode15sMatlabinordinarydifferentialequationnumericalsolutionofapplicationYangHuaZhangLei【Abstract】Manypracticalproblemscanbeusingdifferentialequationsintheformofrepresentation,thetraditionalmethodofsolvingdifferentialequationsaresimilaranalysismethod,tablemethodandgraphicalmethod,thesemethodstocarryonthelargeamountsofhypothesis,sothatthemathematicalmodelhascertaindistortion,havecertainlimitation.Numericalsolutionofusingacomputer,makesolvingmoreaccurateandmoreefficient,andMATLABisakindofmathematicssoftwarepackage,withadvancedprogrammingformat,makingcalculationresultismorecredibility,thereforedifferentialequationandsolutionoftheMATLABinoneoftheapplicationofpracticalsignificance.Thispapernumericalsolutionofdifferentialequationproblemfurtherdiscussion,andtheapplicationofMATLABinwhichthedifficultsolutionimprovementEulermethodandRunge-Kuttamethodontheprogramming,theprogramisconcise,intuitiveandsolutionspeed,methodofpracticalstronger.【Keywords】ordinarydifferentialequation，numericalsolution，Matlab，Eulermethod，Runge-Kuttamethod【引言】微分方程的概念：未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为微分方程是数学科学联系实际问题的主要桥梁之一，它是含有未知数及其导数的方程。常微分方程的求解是现代科学研究和工程技术中经常遇到的问题，然而，从实际问题中建立起来的微分方程往往具有非常复杂的形式，有写解析式难以计算，有的则根本不能用解析式来表达。在实际上对初值问题，一般是要求得到解在若干个点上满足规定精确度的近似值，或者得到一个满足精确度要求的便于计算的表达式，所以利用数值解法求解实际问题显得非常重要。一阶常微分方程的初值，其一般形式为{dydx=f（x，y）y(a)=y0（a≤x≤b）（1）在本文的讨论中，总假设函数f（x，y）连续，且关于y满足莱布尼兹条件，即存在常数L，使得│f(x，y)−f（x，y̅）│≤L|y−y̅|由微分方程理论可知，初值问题式（1）的解必定存在唯一。若给出节点xn=a+nh（n=0,1,2…），其中h为步长，设y（xn）代表方程（1）的精确解在xn的值，yn代表某种算法（忽略计算的舍入误差）得到的y（xn）近似值，所谓数值解法，就是寻求y（x）在系列离散点a=nxxxx210=b处的近似值yn（n=1,2,…），求解过程是顺着节点的次序一步步的向前推进，即按递推公式由已知的y1，y2，y3…yi求出yi+1。1建立数值解法的一些途径1.1、用差商代替导数若步长较小，则有hxyhxyxy)()()('固有公式：0,1,2,n)(),(001xyyyxhfyynnnn此即Euler法。1.2、使用数值积分0),,,',,()(nyyyytF001i)y(xy)f(x,y',1,2,1,0,xynihxi解微分方程：可用以下离散化方法求设对微分方程两边积分得故有公式：)()],(),([200111xyyyxfyxfhyynnnnnn和{yn+1̅̅̅̅̅̅=yn+h(yn−2xn/yn)yn+1=yn+h[(yn−2xnyn)+(yn+1̅̅̅̅̅̅̅−2xnyn+1̅̅̅̅̅̅̅̅)]2此即改进的Euler法。1.3．四阶Runge-Kutta法以Euler法为基础，继续推广和处理，可得一、二、三、四阶Runge-Kutta格式，最常用的一种经典Runge-Kutta格式的具体形式如下：{yn+1=yn+(k1+2k2+2k3+k4)/6k1=hf(xn,yn)k2=hf(xn+h2,yn+k12)k3=hf(xn+h2,yn+k22)k4=hf(xn+h,yn+k3)(n=0,1,2,…)2利用matlab求解一阶常微分方程的初值解在MATLAB中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下：X=dsolve(‘f1’,’f2’,…)函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB具有丰富的函数，我们将其统称为solver，其一般格式为：[T,Y]=solver(odefx,t,y0)该函数表示在区间t=[t0,tf]上，用初始条件y0求解显式常微分方程。solver为命令ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb之一，这些命令各有特点。列表说明如下：求解器ODE类型特点说明ode45非刚性一步算法，4,5阶Runge-Kutta方法累积截断误差大部分场合的首选算法ode23非刚性一步算法，2,3阶Runge-Kutta方法累积截断误差使用于精度较低的情形ode113非刚性多步法，Adams算法，高低精度均计算时间比ode45短))](,())(,([2))(,()()(11111nnnnnnxxnnxyxfxyxfxxdxxyxfxyxynn'(,)yfty3()x3()x可达到ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法，Gear’s反向数值积分，精度中等若ode45失效时，可尝试使用ode23s刚性一步法，2阶Rosebrock算法，低精度。当精度较低时，计算时间比ode15s短odefx为显式常微分方程中的，t为求解区间，要获得问题在其他指定点上的解，则令t=[t0,t1,t2,…]（要求单调），y0初始条件。3根据数值解法的思想编程实现MATLAB中有几个专门用于求解常微分方程的函数,它们的设计思想基于Runge-Kutta方法,基本设计思想为：从改进的欧拉方法比欧拉方法精度高的缘由着手,如果在区间[x1,xi+1]多取几个点的斜率值,然后求取平均值,则可以构造出精度更高的计算方法。这些函数主要包括：ode45、ode23、ode15s、ode113、ode23s、ode23t、ode23tb.其中最常用的是函数ode45,该函数采用变步长四阶五阶Runge-Kutta法求数值解,并采用自适应变步长的求解方法。ode23采用二阶三阶Runge-Kutta法求数值解,与ode45类似,只是精度低一些。ode15s用来求刚性方程组。ode45函数的调用格式为:[tout,yout]=ode45(’yprime’,[t0,tf],y0)其中yprime是表示f(t,y)的M文件名,t0表示自变量的初始值,tf表示自变量的终值,y0表示初始向量值.输出向量tout表示节点(t0;t1;⋯;tn),输出矩阵yout表示数值。例1Rossler方程的数值解法Rossler方程:xzzyxzydtdzdtdydtdx7.52.02.0求解过程如下:(1)建立自定义函数用edit命令建立自定义函数名为rossler.m,内容为:functiondx=rossler(t,x)dx=[-x(2)-x(3);x(1)+0.2*x(2);0.2+(x(1)-5.7)*x(3)];(2)调用对微分方程数值解ode45函数求解用edit命令建立一个命令文件rossler1.m,内容为：x0=[0;0;0];[t,y]=ode45(’rossler’,[0,100],x0);plot(t,y);figure;plot3(y(:,1),y(:,2),y(:,3))得3610~10'(,)yfty(,)fty012,,,tttit如果微分方程由一个或多个高阶常微分方程给出，要得到该方程的数值解,可以将方程转换成一阶常微分方程组。假设高阶常微分方程的一般形式为y(n)=f(t,y,y′,⋯,y(n-1))，而且函数y(t)的各阶导数初值为y(0)，y′(0),…,y(n-1)(0)可以选择一组变量令:x1=y,x2=y′,…,xn=y(n-1),我们就可以把原高阶常微分方程转换成下面的一阶常微分方程组形式：),…,,,(...xxx213221nnxxxtfxx而且初值x1(0)=y(0),x2(0)=y′(0),…，xn(0)=)1(ny(0)。转换以后就可以求原高阶常微分方程的数值解了。例2求微分方程tyeyytytyy22)3(，2)0(y,)0()0(yy的数值解。对方程进行变换,选择变量131212233221321,,txextxxxtxxxxxyxyxyx(1)建立自定义函数用edit命令建立自定义函数名为f.m,内容为:functiony=f(t,x)y=[x(2);x(3);-t^2*x(2)*x(1)^2-t*x(1)*x(3)+exp(-t*x(1))];(2)调用对微分方程数值解ode45函数求解用edit命令建立一个命令文件f2.m,内容为:x0=[2;0;0];[t,y]=ode45(’f’,[0,10],x0);plot(t,y);figure;plot3(y(:,1),y(:,2),y(:,3))得0102030405060708090100-15-10-50510152025-10-5051015-20-100100510152025求刚性微分方程tyeyytytyy22)3(,0)0()0(,2)0(yyy的数值解。调用对微分方程数值解ode15s函数求解x0=[2;0;0];[t,y]=ode15s(’f’,[0,10],x0);plot(t,y(:,1