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

最新推荐文章于 2024-05-14 23:18:50 发布

一夕ξ

最新推荐文章于 2024-05-14 23:18:50 发布

阅读量3.5k

点赞数 3

分类专栏： matlab 文章标签： matlab

本文链接：https://blog.csdn.net/qq_45387575/article/details/119331772

版权

matlab 专栏收录该内容

16 篇文章 1 订阅

订阅专栏

用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可以解出一阶微分方程，或者二阶微分方程中状态变量的解，更加使用与线性系统理论中的状态空间表达式的各个状态变量随时间的变化情况。

一夕ξ

关注

3
点赞
踩
14

收藏

觉得还不错? 一键收藏
打赏
0
评论
matlab--常微分方程的数值解(ODE-s)

用ODE(ordinary differential equations)23和ODE45求解微分方程ODE(ordinary differential equations)23通过使用二阶和三阶隆哥-库塔方法积分微分方程得出
复制链接

扫一扫