matlab--常微分方程的数值解(ODE-s)

用ODE(ordinary differential equations)23和ODE45求解微分方程

ODE(ordinary differential equations)23

[t,y]=ode23('func_name',[start_time,end_time],y(0))

通过使用二阶和三阶隆哥-库塔方法积分微分方程得出

function ydot=eq1(t,y)
ydot=cos(t)

[t,y]=ode23('eq1',[0 2*pi],2);----------就会产生一系列的数值解

ydot =

     1


ydot =

    0.9968


ydot =

    0.9928


ydot =

    0.9872


ydot =

    0.9400


ydot =

    0.9038


ydot =

    0.8596


ydot =

    0.7394


ydot =

    0.6676


ydot =

    0.5890


ydot =

    0.3909


ydot =

    0.2835


ydot =

    0.1724


ydot =

   -0.1176


ydot =

   -0.2604


ydot =

   -0.3977


ydot =

   -0.6618


ydot =

   -0.7709


ydot =

   -0.8610


ydot =

   -0.6018


ydot =

   -0.6919


ydot =

   -0.7723


ydot =

   -0.8998


ydot =

   -0.9450


ydot =

   -0.9770


ydot =

   -1.0000


ydot =

   -0.9953


ydot =

   -0.9799


ydot =

   -0.9248


ydot =

   -0.8846


ydot =

   -0.8365


ydot =

   -0.7234


ydot =

   -0.6577


ydot =

   -0.5865


ydot =

   -0.4280


ydot =

   -0.3429


ydot =

   -0.2550


ydot =

   -0.0510


ydot =

    0.0524


ydot =

    0.1553


ydot =

    0.4529


ydot =

    0.5868


ydot =

    0.7063


ydot =

    0.3516


ydot =

    0.4448


ydot =

    0.5334


ydot =

    0.6932


ydot =

    0.7628


ydot =

    0.8245


ydot =

    0.9143


ydot =

    0.9477


ydot =

    0.9730


ydot =

    0.9932


ydot =

    0.9983


ydot =

     1

>> [t,y]=ode23('eq1',[0 2*pi],2);

ydot =

     1


ydot =

    0.9968


ydot =

    0.9928


ydot =

    0.9872


ydot =

    0.9400


ydot =

    0.9038


ydot =

    0.8596


ydot =

    0.7394


ydot =

    0.6676


ydot =

    0.5890


ydot =

    0.3909


ydot =

    0.2835


ydot =

    0.1724


ydot =

   -0.1176


ydot =

   -0.2604


ydot =

   -0.3977


ydot =

   -0.6618


ydot =

   -0.7709


ydot =

   -0.8610


ydot =

   -0.6018


ydot =

   -0.6919


ydot =

   -0.7723


ydot =

   -0.8998


ydot =

   -0.9450


ydot =

   -0.9770


ydot =

   -1.0000


ydot =

   -0.9953


ydot =

   -0.9799


ydot =

   -0.9248


ydot =

   -0.8846


ydot =

   -0.8365


ydot =

   -0.7234


ydot =

   -0.6577


ydot =

   -0.5865


ydot =

   -0.4280


ydot =

   -0.3429


ydot =

   -0.2550


ydot =

   -0.0510


ydot =

    0.0524


ydot =

    0.1553


ydot =

    0.4529


ydot =

    0.5868


ydot =

    0.7063


ydot =

    0.3516


ydot =

    0.4448


ydot =

    0.5334


ydot =

    0.6932


ydot =

    0.7628


ydot =

    0.8245


ydot =

    0.9143


ydot =

    0.9477


ydot =

    0.9730


ydot =

    0.9932


ydot =

    0.9983


ydot =

     1

>> f=2+sin(t);
>> plot(t,y,'o',t,f)

这一系列数值解就是原函数的曲线上。y是数值解,f是解析解

查看数值解与解析解的误差

 >> for i=1:1:size(y)
err(i)=abs((f(i)-y(i))/f(i));
end

ODE(ordinary differential equations)45

[t,w]=ode45('function',[start_time,end_time],y(0))

使用高阶龙格-库塔公式

>> [t,w]=ode45('eq1',[0 2*pi],2);

ydot =

     1


ydot =

    0.9968


ydot =

    0.9927


ydot =

    0.9488


ydot =

    0.9369


ydot =

    0.9203


ydot =

    0.9203


ydot =

    0.8640


ydot =

    0.8307


ydot =

    0.6180


ydot =

    0.5732


ydot =

    0.5146


ydot =

    0.5146


ydot =

    0.4031


ydot =

    0.3449


ydot =

    0.0379


ydot =

   -0.0179


ydot =

   -0.0876


ydot =

   -0.0876


ydot =

   -0.2118


ydot =

   -0.2727


ydot =

   -0.5567


ydot =

   -0.6022


ydot =

   -0.6564


ydot =

   -0.6564


ydot =

   -0.7458


ydot =

   -0.7862


ydot =

   -0.9387


ydot =

   -0.9564


ydot =

   -0.9745


ydot =

   -0.9745


ydot =

   -0.9949


ydot =

   -0.9993


ydot =

   -0.9621


ydot =

   -0.9454


ydot =

   -0.9203


ydot =

   -0.9203


ydot =

   -0.8640


ydot =

   -0.8307


ydot =

   -0.6180


ydot =

   -0.5732


ydot =

   -0.5146


ydot =

   -0.5146


ydot =

   -0.4031


ydot =

   -0.3449


ydot =

   -0.0379


ydot =

    0.0179


ydot =

    0.0876


ydot =

    0.0876


ydot =

    0.2118


ydot =

    0.2727


ydot =

    0.5567


ydot =

    0.6022


ydot =

    0.6564


ydot =

    0.6564


ydot =

    0.7458


ydot =

    0.7862


ydot =

    0.9387


ydot =

    0.9564


ydot =

    0.9745


ydot =

    0.9745


ydot =

    0.9836


ydot =

    0.9875


ydot =

    0.9990


ydot =

    0.9997


ydot =

     1


ydot =

     1

>> f=2+sin(t);
>> plot(t,w,'0',t,f)

 


ode45返回了45个数据点,而ode23返回7个数据点。


求两个一阶方程组的数值解:

\frac{\mathrm{dx} }{\mathrm{d} t}=-x^2+y

\frac{\mathrm{dy} }{\mathrm{d} t}=-x+xy

.m文件:

function xdot=eqx(t,x)
xdot=zeros(2,1)
xdot(1)=-x(1)^2+x(2);
xdot(2)=-x(1)-x(1)*x(2)

>> [t,x]=ode45('eqx',[0 10],[0 1]);

>> plot(t,x(:,1),t,x(:,2),'--')

 >> plot(x(:,1),x(:,2)),xlabel('x(1)'),ylabel('x(2)')

求解二阶微分方程

{y}''+16y=sin(4.3t)

y(0)=y'(0)=0


首先转换成一个一阶微分方程组:

x_1=y,x_2=y'

x_1'=y'=x_2;x_2'=y''=sin(4.3t)-16x_1


.m文件

function xdot=eqx(t,x)
xdot=zeros(2,1)
xdot(1)=x(2);
xdot(2)=sin(4.3*t)-16*x(1);

命令行窗口:

>> [t,x]=ode45('eqx',[0 2*pi],[0,0]);

>> plot(t,x(:,1),t,x(:,2),'--')

 >> plot(x(:,1),x(:,2)),xlabel('x(1)'),ylabel('x(2)')

 


求解:

x_1'=2x_1x_2;x_2'=-x_1^2;x_1(0)=x_2(0)=0


.m文件:

function xdot=eqx(t,x)
xdot=zeros(2,1)----------------------------(分配存储空间)
xdot(1)=2*x(1)*x(2);
xdot(2)=x(1)^2;

>>[t,x]=ode45('eqx',[0 2*pi],[0,0])

>> plot(t,x(:,1),t,x(:,2),'--')

 

>>[t,x]=ode45('eqx',[0 2*pi],[1,2])

>> plot(t,x(:,1),t,x(:,2),'--')

 


求解:

y'=-2.3y;y(0)=0


function ydot=eqx(t,y)
ydot=-2.3*y;

>> [t,x]=ode45('eqx',[0 2*pi],[1,2])

>> plot(t,y)

 >> syms y(t);
>> eqn=diff(y,t)+2.3*y==0;
>> cond=(y(0)==0);
>> s=dsolve(eqn,cond)
 
s =
 
0


求解:

y''-2y'+y=exp(-t);y(0)=2,y'(0)=0



x_1=y;x_2=y'

x_1'=y'=x_2

x_2'=y''=exp(-t)-y+2y'=exp(-t)-x_1+2x_2


.m文件:

>> [t,x]=ode45('eqx',[0 2*pi],[2,0])

>> plot(t,x(:,1),t,x(:,2),'--')

 >> syms y(t);
>> eqn=diff(y,t,2)-2*diff(y,t)+y==exp(-t);
>> Dy=diff(y,t);
>> cond=[y(0)==2,Dy(0)==0];
>> s=dsolve(eqn,cond)
 
s =
 
(exp(-t)*(7*exp(2*t) - 6*t*exp(2*t) + 1))/4

>> ezplot(s)

 


直接用dsolve:可以直接解出方程

而用ode45可以解出一阶微分方程,或者二阶微分方程中状态变量的解,更加使用与线性系统理论中的状态空间表达式的各个状态变量随时间的变化情况。


 

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

一夕ξ

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值