第3-1章 连续系统数值积分法仿真Matlab编程

% 输出参数：NewX是这一步计算的新的状态向量
Fra Baidu bibliotek
单步数值积分函数只是对微分方程组StateModel进行一步的计算， 计算法由InteMethod参数指定，可以上欧拉法，RK2或RK4。
4
function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel); if InteMethod == 'RK4' k1=eval([StateModel,'(curt,curx,curu)']); k2=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k1,curu)']); k3=eval([StateModel,'(curt+0.5*h,curx+0.5*h*k2,curu)']); k4=eval([StateModel,'(curt+h,curx+h*k3,curu)']); NewX=curx+h*(k1+2*k2+2*k3+k4)/6; elseif InteMethod == 'RK2' k1=eval([StateModel,'(curt,curx,curu)']); k2=eval([StateModel,'(curt+h,curx+h*k1,curu)']); NewX=curx+0.5*h*(k1+k2); else %欧拉法EUL Qk=eval([StateModel,'(curt,curx,curu)']);%eval用来执行一个函数，传递函数 名和函数的输入参数 NewX=curx+h*Qk; end
6
3.2 算法稳定性分析仿真编程
针对下面的系统，求用解析法、欧拉法和RK4分别求解，计算欧 拉法最大允许步长，将步长从0.1逐渐增大，比较三种解的效果。
4 x x x (0) 1
解：1）该系统是稳定的，解析解为
(3.2-1)
x(t) e4t
2）用欧拉法计算本例时，其步长应该满足
第三章 Matlab编程(数值积分法仿真)
3.1 连续系统数值积分法仿真编程思路
目的：针对下面的系统编写通用的数值积分法仿真程序
u g (t , x , u , y ) f (t , x , u ) x y h (t , x , u )
给定初值：u0，y0,系统中的向量都采用列向量表示
5
函数w_DigiInteSimu和w_StepIntegral构造了一个数值积分法 仿真的框架，并不涉及具体的系统。
具体的系统由StateModel,ControlFile,OutputFile参个参数决定， 实际上就是三个函数文件名，这三个函数输入输出参数必须遵循特定 的格式，在准备好由这3个函数描述的系统后，调用w_DigiInteSimu 即可进行仿真。 还需要一个调用w_DigiInteSimu进行仿真的程序，它指定模型文 件，指定初始参数，并且对仿真结果绘图。
(3.1-1)
对这样的系统进行仿真，实际上涉及到控制的计算、状态的数值 积分计算和输出的计算。3个函数g,f,h确定后，就可以完整地描 述一个系统。
1
(1)数值积分法仿真框架函数
function [t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile, ControlFile) % 函数功能：用数值积分法仿真 % 输入参数：tstart, tstop,h 分别是起始时间、结束时间和仿真步长，是标量 % x0,u0是状态、输入的初值，都是列向量 % cnty是输出变量的个数 % InteMethod时数值积分方法，可以是'EUL' ,'RK2','RK4' % StateModel是状态模型的文件名 % ControlFile是控制作用的文件名 % OutputFile是系统输出的文件名 % 输出参数：t是仿真结果的时间序列 % y是仿真结果系统的输出序列
2
function [t,y]=w_DigiInteSimu(tstart,tstop,h,x0,u0,cnty,InteMethod,StateModel,OutputFile, ControlFile) t=[tstart:h:tstop];%t数一个行序列 cntt=size(t,2);%返回列数 y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果 y0=eval([OutputFile,'(tstart,x0,u0)']);%计算初始输出 y(:,1)=y0’;%将cury作为输出的第1列 curx=x0; %当前一步的x curu=u0; %当前一步的u cury=y0; %当前一步的y for i=1:1:cntt-1 curu=eval([ControlFile,'(t(i),h,curx,curu,cury)']);%计算控制时传递的参数：当前 时间，步长，当前状态和输出 curx=w_StepIntegral(t(i),h,curx,curu,InteMethod,StateModel);%单步积分运算 cury=eval([OutputFile,'(t(i),curx,curu)']);%计算输出 y(:,i+1)=cury‘;%将输出加入到输出序列里 end
1 4 h 1 h 0 . 5
3）RK4步长条件式是一个高阶不等式，无法直接求解，只能用 试探法确定RK4的步长上限。
3
(2)单步数值积分法函数
function NewX=w_StepIntegral(curt,h,curx,curu,InteMethod,StateModel) %函数功能：单步积分运算,与模型方程有关 % 输入参数：curt是当前时间，h是数值积分步长 % % % curx,curu分别是当前的状态和控制向量 InteMethod是积分算法，字符串类型，可以取'EUL','RK2','RK4' StateModel是状态模型文件名称，字符串类型