用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个数据点。
求两个一阶方程组的数值解：
.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)')
求解二阶微分方程
首先转换成一个一阶微分方程组：
.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)')
求解：
.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),'--')
求解：
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
求解：
.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可以解出一阶微分方程，或者二阶微分方程中状态变量的解，更加使用与线性系统理论中的状态空间表达式的各个状态变量随时间的变化情况。
