数值计算(1) --求解连续微分系统和混沌系统

前言

微分系统在工程项目中很常见，通过物理建模之后，基本都需要求解微分方程得到其结果，混沌系统属于特殊的一类微分系统，在某些项目上也很常见，同时可以引申出分岔图、李雅普诺夫指数谱、相图、庞加莱截面等，本文探讨通过matlab常见的微分求解函数和simulink求解器来实现计算。

关键字：微分系统，混沌系统，Simulink

正文

1、常微分方程（Lorenze混沌系统）

方法1：m文件实现

x0=[0;0;1e-3]; %设定初始值[t,x]=ode45(@lorenzfun,[0,100],x0); %调用函数ode45求解，figure(1)plot(t,x)figure(2)plot3(x(:,1),x(:,2),x(:,3))

function dx=lorenzfun(t,x)% 输入微分方程a=10;c=28;b=8/3;dx=zeros(3,1);dx(1)=-b*x(1)+x(2)*x(3);dx(2)=-a*x(2)+10*x(3);dx(3)=-x(1)*x(2)+c*x(2)-x(3);

结果如图

方法2：Simulink模块实现

其中三个积分模块的初始值设置与exam1相同，仿真时长为100s。精度设置：Simulation--Configuration Parameters—Relative tolerance, 1e-3改为1e-5（试试不作此修改的结果比较）。运行后双击示波器scope后可看到。

在matlab命令窗口输入画图命令：

figureplot(tout,yout)figureplot3(yout(:,2),yout(:,3),yout(:,1))

方法3：simulink向量模块

在Fcn模块里面分别定义好3组微分方程，最后进行积分求解即可

2、常时滞微分方程

方法1：m文件需调用dde23来求解

sol = dde23('exam1f',[1, 0.2],ones(3,1),[0, 5]);plot(sol.x,sol.y);title('Example 2')xlabel('time t');ylabel('y(t)');

function v = exam1f(t,y,Z)ylag1 = Z(:,1);ylag2 = Z(:,2);v = zeros(3,1);v(1) = ylag1(1);v(2) = ylag1(1) + ylag2(2);v(3) = y(2);

方法2：Simulink中S函数来实现

注：用Simulink中S函数求解时滞微分方程的核心思想在于：将时滞变量作为S函数的外部输入，这个需要通过transport delay模块实现。

延申思考

1、在求解微分方程后如何得到分叉图？

Tips:系统单参数分岔图的计算方法：最大值法和Poincare截面法，最大值法最为简便，对系统微分方程（组）进行求解，对求解的结果用getmax函数进行取点，并绘图即可。