Matlab通过ode系列函数求解微分方程

✅作者简介：热爱科研的算法开发者，Python、Matlab项目可交流、沟通、学习。

🍎个人主页：算法工程师的学习日志

MATLAB有很多用于求解微分方程的内置函数。MATLAB包含了用于求解常微分方程（ODE）的函数，微分表达式一般如下

对于高阶微分方程必须重新表述为一个一阶系统微分方程。

并不是所有的微分方程都可以用同样的方法求解，所以MATLAB提供了许多不同的常微分方程求解器，如ode45、ode23、ode113等。

考虑一个细菌种群数学模型，x为现在的细菌数量，细菌生长率为bx，死亡率为px^2，其数学表达式为：

其中 b=1，p=0.5

登录后复制

         
         

          
          function dx = bacteriadiff(t,x)
          
          
b=1;
          
          
p=0.5;
          
          
dx = b*x - p*x^2;
         
         

• 1.
• 2.
• 3.
• 4.

登录后复制

         
         

          
          clear
          
          
clc
          
          
tspan=[0 1];
          
          
x0=100;
          
          
[t,y]=ode45(@bacteriadiff, tspan,x0);
          
          
plot(t,y)
         
         

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.

对于采用变参数的微分数学模型方法

 其中，假定a = 1/T,T是仿真的时间，b = 1,x(0) = 1, T = 5

 

登录后复制

         
         

          
          function dx = mysimplediff(t,x,param)
          
          
a=param(1);
          
          
b=param(2);
          
          
dx=a*x+b;
         
         

• 1.
• 2.
• 3.
• 4.

登录后复制

         
         

          
          tspan=[0 25];
          
          
x0=1;
          
          
a=-1/5;
          
          
b=1;
          
          
param=[a b];
          
          
[t,y]=ode45(@mysimplediff, tspan, x0,[], param);
          
          
plot(t,y)
         
         

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.

使用ode23函数求解微分方程并绘制[t0,tf]区间上

 假定

微分方程可表达为：

登录后复制

         
         

          
          function dw = diff_task3(t,w)
          
          
dw = -(1.2 + sin(10*t))*w;
         
         

• 1.
• 2.

tspan=[0 5];

w0=1;

[t,w]=ode23(@diff_task3, tspan, w0);

plot(t,w)

 

求解含有二阶的微分方程

令：

高阶的系统（二阶、三阶等）需要降为一阶来书写表达式，学过现代控制理论的应该熟悉这个

令：

则

登录后复制

         
         

          
          function dx = diff_secondorder(t,x) 
          
          
[m,n] = size(x);
          
          
dx = zeros(m,n)
          
          
dx(1) = x(2);
          
          
dx(2) = (2-2*t*x(2)-3*x(1))/(1+t^2);
         
         

• 1.
• 2.
• 3.
• 4.
• 5.

登录后复制

         
         

          
          tspan=[0 5];
          
          
x0=[0; 1];
          
          
[t,x]=ode23(@diff_secondorder, tspan, x0);
          
          
plot(t,x)
          
          
legend('x1','x2')
          
          
tspan=[0 5];
          
          
x0=[0; 1];
          
          
[t,x]=ode23(@diff_secondorder, tspan, x0);
          
          
plot(t, x(:,2))
         
         

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.
• 8.
• 9.

原创作者: domi 转载于: https://blog.51cto.com/domi/6139313