资源描述:
《matlab求解常微分方程new》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、用matlab求解常微分方程 在MATLAB中,由函数dsolve()解决常微分方程(组)的求解问题,其具体格式如下:r=dsolve('eq1,eq2,...','cond1,cond2,...','v')'eq1,eq2,...'为微分方程或微分方程组,'cond1,cond2,...',是初始条件或边界条件,'v'是独立变量,默认的独立变量是't'。函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解。例1:求解常微分方程的MATLAB程序为:dsolve('Dy=1/(x+y)','x'
2、),注意,系统缺省的自变量为t,因此这里要把自变量写明。其中:Y=lambertw(X)表示函数关系Y*exp(Y)=X。例2:求解常微分方程的MATLAB程序为:Y2=dsolve('y*D2y-Dy^2=0','x')Y2=dsolve('D2y*y-Dy^2=0','x')我们看到有两个解,其中一个是常数0。例3:求常微分方程组通解的MATLAB程序为:[X,Y]=dsolve('Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)','t')例4:求常微分方程组通解的MATLAB程序为:[X,Y]=dsolve('Dx+2
3、*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t)','x(0)=2,y(0)=0','t')以上这些都是常微分方程的精确解法,也称为常微分方程的符号解。但是,我们知道,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB具有丰富的函数,我们将其统称为solver,其一般格式为:[T,Y]=solver(odefun,tspan,y0)该函数表示在区间tspan=[t0,tf]上,用初始条件y0求解显式常微分方程。solver为命
4、令ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一,这些命令各有特点。我们列表说明如下:求解器特点说明ode45一步算法,4,5阶Runge-Kutta方法累积截断误差大部分场合的首选算法ode23一步算法,2,3阶Runge-Kutta方法累积截断误差使用于精度较低的情形ode113多步法,Adams算法,高低精度均可达到计算时间比ode45短ode23t采用梯形算法适度刚性情形ode15s多步法,Gear’s反向数值积分,精度中等若ode45失效时,可尝试使用ode23s一步法,2阶Rose
5、brock算法,低精度。当精度较低时,计算时间比ode15s短odefun为显式常微分方程中的tspan为求解区间,要获得问题在其他指定点上的解,则令(要求单调递增或递减),y0初始条件。例5:求解常微分方程,,的MATLAB程序如下:y=dsolve('Dy=-2*y+2*x^2+2*x','y(0)=1','x')x=0:0.01:0.5;yy=subs(y,x);fun=inline('-2*y+2*x*x+2*x');[x,y]=ode15s(fun,[0:0.01:0.5],1);ys=x.*x+exp(-2*x);plot(x,y,
6、'r',x,ys,'b')例6:求解常微分方程的解,并画出解的图形。分析:这是一个二阶非线性方程(函数以及所有偏导数军委一次幂的是现性方程,高于一次的为非线性方程),用现成的方法均不能求解,但我们可以通过下面的变换,将二阶方程化为一阶方程组,即可求解。令:,,,则得到:解:function[dfy]=mytt(t,fy)%f1=y;f2=dy/dt%求二阶非线性微分方程时,把一阶、二阶直到(n-1)阶导数用另外一个函数代替%用ode45命令时,必须表示成Y'=f(t,Y)的形式%Y=[y1;y2;y3],Y'=[y1';y2';y3']=[y2
7、;y3;f(y1,y2,y3)],%其中y1=y,y2=y',y3=y''%更高阶时类似dfy=[fy(2);7*(1-fy(1)^2)*fy(2)-fy(1)];clear;clc[t,yy]=ode45('mytt',[040],[1;0]);plot(t,yy)legend('y','dy') 【例4.14.2.1-1】采用ODE解算指令研究围绕地球旋转的卫星轨道。(1)问题的形成轨道上的卫星,在牛顿第二定律,和万有引力定律作用下有,引力常数G=6.672*10-11(N.m2/kg2),ME=5.97*1024(kg)是地球的质量。假定
8、卫星以初速度vy(0)=4000m/s在x(0)=-4.2*107(m)处进入轨道。(2)构成一阶微分方程组令Y=[y1y2y3y4]T=[xyvxv