[讨论]

multism

wangrongdong88

2014-4-1 19:22:48 4178

用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') ，

注意，系统缺省的自变量为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*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为命令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阶Rosebrock算法，低精度。	当精度较低时，计算时间比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,'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;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',[0 40],[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)是地球的质量。假定卫星以初速度vy(0)=4000m/s在x(0)=-4.2*107(m)处进入轨道。

（2）构成一阶微分方程组

令Y=[y1 y2y3 y4]T=[x y vx vy]T=[x y x' y']T

（3）根据上式

[dYdt.m]

function Yd=DYdt(t,Y)

% t

% Y

global GME %

xy=Y(1:2);Vxy=Y(3:4); %

r=sqrt(sum(xy.^2));

Yd=[Vxy;-G*ME*xy/r^3]; %

（4）

global G ME % <1>

G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;

tspan=[t0,tf]; %

Y0=[x0;0;0;vy0]; %

[t,YY]=ode45('DYdt',tspan,Y0);% <8>

X=YY(:,1); %

Y=YY(:,2); %

plot(X,Y,'b','Linewidth',2); hold on

%axis('image') %

[XE,YE,ZE] = sphere(10); %

RE=0.64e7; %

XE=RE*XE;YE=RE*YE;ZE=0*ZE; %

mesh(XE,YE,ZE),hold off %

练习：

1．利用MATLAB求常微分方程的初值问题，的解。

r=dsolve('Dy+3*y=0','y(0)=2','x')

2．利用MATLAB求常微分方程的初值问题，，的解。

r=dsolve('D2y*(1+x^2)-2*x*Dy=0','y(0)=1,Dy(0)=3','x')

3．利用MATLAB求常微分方程的解。

解：y=dsolve('D4y-2*D3y+D2y','x')

4．利用MATLAB求常微分方程组的特解。

[X,Y]=dsolve('Dx*2+4*x+Dy-y=exp(t),Dx+3*x+y=0','x(0)=1.5,y(0)=0','t')

5．求解常微分方程，，，的特解，并作出解函数的曲线图。

r=dsolve('D2y-2*(1-y^2)*Dy+y=0','y(0)=0,Dy(0)=0','x')

function DY=mytt2(t,Y)

DY=[Y(2);2*(1-Y(1)^2)*Y(2)+Y(1)];

clear;clc

[t,YY]=ode45('mytt2',[030],[1;0]);

plot(t,YY)

legend('y','dy')

求解特殊函数方程

勒让德方程的求解

r=dsolve('(1-x^2)*D2y-2*x*Dy+y*l*(l+1)=0','x')

连带勒让德方程的求解：

r=dsolve('(1-x^2)*D2y-2*x*Dy+y*(l*(l+1)-m^2/(1-x^2))=0','x')

贝塞尔方程

r=dsolve('x^2*D2y+x*Dy+(x^2-v^2)*y=0','x')

利用maple

mapledsolve((1-x^2)*diff(y(x),x$2)-2*x*diff(y(x),x)+n*(n+1)*y(x) = 0, y(x))

利用MAPLE的深层符号计算资源经典特殊函数的调用MAPLE库函数在线帮助的检索树发挥MAPLE的计算潜力调用MAPLE函数

【例5.7.3.1-1】求递推方程的通解。

（1）

gs1=maple('rsolve(f(n)=-3*f(n-1)-2*f(n-2),f(k));')

（2）调用格式二

gs2=maple('rsolve','f(n)=-3*f(n-1)-2*f(n-2)','f(k)')

【例5.7.3.1-2】求的Hessian矩阵。

（1）

FH1=maple('hessian(x*y*z,[x,y,z]);')

（2）

FH2=maple('hessian','x*y*z','[x,y,z]')

（3）

FH=sym(FH2)

【例5.7.3.1-3】求在处展开的截断8阶小量的台劳近似式。

（1）

TL1=maple('mtaylor(sin(x^2+y^2),[x=0,y=0],8)')

（2）

maple('readlib(mtaylor);');

TL2=maple('mtaylor(sin(x^2+y^2),[x=0,y=0],8)');

pretty(sym(TL2))

运行MAPLE程序

【例5.7.3.2-1】目标：设计求取一般隐函数的导数解析解的程序，并要求：该程序能象MAPLE原有函数一样可被永久调用。

（1）

[DYDZZY.src]

DYDZZY:=proc(f)

# DYDZZY(f) is used to get the derivateof

# an implicit function

local Eq,deq,imderiv;

Eq:='Eq';

Eq:=f;

deq:=diff(Eq,x);

readlib(isolate);

imderiv:=isolate(deq,diff(y(x),x));

end;

（2）

procread('DYDZZY.src')

（3）

s1=maple('DYDZZY(x=log(x+y(x)));')

s2=maple('DYDZZY(x^2*y(x)-exp(2*x)=sin(y(x)))')

s3=maple('DYDZZY','cos(x+sin(y(x)))=sin(y(x))')

（4）

clear maplemex

procread('DYDZZY.src');

maple('save(`DYDZZY.m`)');

（5）

maple('read','`DYDZZY.m`');

ss2=maple('DYDZZY(x^2*y(x)-exp(2*x)=sin(y(x)))')

2014-4-1 19:22:48　　评论 淘帖0 举报

撰写讨论

你正在撰写讨论

如果你是对讨论或其他讨论精选点评或询问，请使用“评论”功能。

声明：本文内容及配图由入驻作者撰写或者入驻合作网站授权转载。文章观点仅代表作者本人，不代表电子发烧友网立场。文章及其配图仅供工程师学习之用，如有内容图片侵权或者其他问题，请联系本站作侵删。侵权投诉

发讨论

精选推荐

multisim14.0，变压器仿真为什么出现这样的错误结果？

11078 浏览 2 评论
NCP1015反激开关电源multisim仿真错误

14944 浏览 1 评论
Multisim14.2中CD4538高电平输出为什么只有5V？

17813 浏览 2 评论
分析电流检测运放电路，输入电流I和输出电压Vout的关系

16736 浏览 1 评论
如何使BT33触发电路输出可调？

24498 浏览 4 评论

快速回复 返回顶部 返回列表

关注微信公众号

电子发烧友网

电子发烧友论坛

社区合作: 刘勇; 联系电话：15994832713; 邮箱地址：liuyong@huaqiu.com

社区管理: elecfans短短; 微信：elecfans_666; 邮箱：users@huaqiu.com

【优惠升级】华秋PCB首单最高立减100元，SMT免费贴片！！！

返回 multisim论坛

搜索历史

multism

相关推荐

撰写讨论

精选推荐

站长推荐 /7