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

最新推荐文章于 2024-09-14 23:09:27 发布

一夕ξ

最新推荐文章于 2024-09-14 23:09:27 发布

阅读量3.6k

点赞数 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)