控制系统数字仿真 四阶龙格库塔法
数值计算课程设计四阶RungeKutta方法
湖南工业大学课程设计资料袋理学院（系、部）2013 学年第 2 学期课程名称数值计算方法指导教师职称副教授学生姓名专业班级信息与计算科学班学号学生姓名专业班级信息与计算科学1002班学号学生姓名专业班级信息与计算科学学号题目四阶Runge-Kutta方法成绩起止日期2013 年6 月24日～2013 年7月5日目录清单湖南工业大学课程设计任务书2012 —2013 学年第2 学期理学院（系、部）信息与计算科学专业1002 班级课程名称：数值计算方法设计题目：四阶Runge-Kutta方法完成期限：自2013 年 6 月24 日至2013 年7月 5 日共 2 周指导教师（签字）：年月日系（教研室）主任（签字）：年月日数值计算方法设计说明书四阶Runge-Kutta方法起止日期：2013 年6 月24 日至2013 年7月 5 日学生姓名班级信息与计算科学班学号成绩指导教师(签字)理学院（院、部）2013年7月5日目录一、摘要 (5)二、问题重述 (5)三、方法原理及实现 (5)四、计算公式或算法 (5)五、Matlab程序 (6)六、测试数据及结果 (6)七、结果分析 (10)八、方法改进 (10)九、心得体会 (10)十、参考文献 (10)一、摘要本课程设计主要内容是用四阶Runge-Kutta 方法解决常微分方程组初值问题的数值解法，通过分析给定题目使用Matlab 编写程序计算结果并绘图，最后对计算结果进行分析，得到结论。
二、问题重述在计算机上实现用四阶Runge-Kutta 求一阶常微分方程初值问题()()[]()⎩⎨⎧=∈=1,,,,'y a y b a x y x f x y的数值解，并利用最后绘制的图形直观分析近似解与准确解之间的比较。
三、方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
matlab经典的4级4阶runge kutta法 -回复
matlab经典的4级4阶runge kutta法-回复使用MATLAB 实现经典的4 阶4 级Runge-Kutta 法引言：数值计算是现代科学和工程中的一个重要领域，它涉及到通过计算机模拟来解决数学问题。
在数值计算中，求解微分方程是一个常见的任务。
Runge-Kutta 法是求解微分方程的一种常见方法，它可以用于数值求解常微分方程和偏微分方程。
本文将介绍经典的4 级4 阶Runge-Kutta 法的原理，并使用MATLAB 来实现该方法。
一、原理介绍：Runge-Kutta 法是数值计算领域中最常用的方法之一。
它通过将微分方程的解逐步逼近来求解微分方程。
经典的4 级4 阶Runge-Kutta 法基于以下公式：\begin{align*}k_1 &= h f(t_n, y_n) \\k_2 &= h f(t_n + \frac{h}{2}, y_n + \frac{k_1}{2}) \\k_3 &= h f(t_n + \frac{h}{2}, y_n + \frac{k_2}{2}) \\k_4 &= h f(t_n + h, y_n + k_3) \\y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)\end{align*}其中，h 是步长，t_n 是当前时间点，y_n 是当前的解，f(t, y) 是微分方程的右手函数。
二、算法实现：现在我们将使用MATLAB 实现经典的4 级4 阶Runge-Kutta 法，并解决一个简单的一阶常微分方程。
首先，我们定义一个MATLAB 函数，用于实现4 级4 阶Runge-Kutta 法。
函数接受输入参数为微分方程的右手函数f(t, y)，初始时间t_0，初始解y_0，以及步长h。
函数输出为一个数组，包含了每个时间点的解。
以下是MATLAB 代码实现：matlabfunction y = runge_kutta(f, t0, y0, h, num_steps)初始化解数组y = zeros(num_steps+1, 1);y(1) = y0;循环计算每个时间点的解for i = 1:num_stepst = t0 + (i-1)*h;计算k1, k2, k3, 和k4k1 = h * f(t, y(i));k2 = h * f(t + h/2, y(i) + k1/2);k3 = h * f(t + h/2, y(i) + k2/2);k4 = h * f(t + h, y(i) + k3);计算下一个时间点的解y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4)/6;endend接下来，我们使用这个函数来解决一个简单的一阶常微分方程。
控制系统数字仿真模拟题
控制系统数字仿真模拟题一、填空题1、数值积分法中，计算精度p=2的是 梯形法2、 混合法 是机理模型法和统计模型法的结合3、柔性制造系统属于 离散事件 系统.4、零极点增益形式可用于分析系统的 稳定性 和 快速性5、 现实性 、 简洁性 、 适应性 是建立系统模型应该依照的原则.6、系统的三大要素为： 实体 、 属性 和活动。
7、通常仿真时多采用四阶龙格 库塔法 、其原因就是这种计算公式的截断误差较小.8、 相似论 是系统仿真的主要依据.9、一个电机转速控制系统中，属于电机所具有的属性的为: 电机转速10、我们在选择数值算法的时候要 考虑精度 、 计算速度 以及稳定性等原则进行.二、单选题1.运行下列命令后A1=[1，2，3；4，5，6；7，8，9]；A2=A1；A3=cat(1，A1，A2)，系统输出结果为( B )A.123147456258789369B.123456789147258369C.123456789D.147258369 2.设某一系统的状态方程矩阵为a=[-3，1；1，-3]；b=[1，1；1，1]；c=[1，1；1，-1]；d=[0]并且执行后得可控性矩阵和可观性矩阵的秩分别为cam=ctrb(a ，b)＝1，rcam=rank(cam)＝2，因此这一系统为( D )A.不可控且不可观的系统B.可控且可观的系统C.可控但不可观的系统D.不可控但可观的系统3.可以将模块按照顺时针进行旋转的快捷键为( A )A.ctrl+rB.ctrl+yC.alt+rD.alt+y4.在Matlab 系统中,调用Simulink环境的工具栏图标为( A )A. B. C.5.下列符号中可以引导注释行的是( D )A.&B.@C.$D.%6.若A=412303214--⎡⎤⎢⎥-⎢⎥⎢⎥-⎣⎦，则C=(A>0)&(A<3)的结果为( B )A.001001011B.001000010C.111110110D.0011100107.MATLAB系统中若要使系统选择short和shortE中最好的表示，则采用命令( C )A.shortB.shortEC.shortGD.longE8.列出工作内存中的变量名称以及细节，只需在命令窗口输入( A )A.whatB.whoC.echoonD.whose9.设一个五阶魔方阵B=magic(5)，提取B阵的第1行，第2行的第1，3，5个元素的命令为( B )A.B(1，2:[1，3，5])B.B([1:2]，[1，3，5])C.B([1:2]，1:3:5))D.B(1:2；[1，3，5])10.下列命令中可以创建起始值为0，增量值为0.5，终止值为10的等差数列的是( A )A.a＝0:0.5:10B.a＝linspace(0，10，0.5)C.linspace(0，10，10)D.logspace(0，1，11)11.若a=[102；300；130；111]，则any(a)=( C )A.011B.110C.111D.10012.设s＝‘haha’，可以看到字符s的ascii码值的命令为( C )A.size(s）B.isstr(s）C.abs(s）D.eval(s）13.PSPICE是( B )软件.A.模型及混合信号仿真软件B.模拟电路仿真软件C.机械系统动力学自动分析软件D.大型通用有限元分析软件14.将多项式2(22)(4)(1)s s s s++++展开的命令中正确的是( D )A.conv([1，2，2]，conv([4，1]，[1，1]))B.conv([2，2，1]，conv([4，1]，[1，1]))C.conv([2，2，1]，conv([1，4]，[1，1]))D.conv([1，2，2]，conv([1，4]，[1，1]))15.w=conv([1，2，3]，conv([1，2]，[1，1]))的值为( C )A.3111372B.2713113C.1511136D.151422219三、判断题1．影响系统而又不受系统直接控制的全部外界因素的集合叫外部活动.( 错）2．系统仿真就是建立系统的动态模型并在模型上进行实验(或试验).（对）3．状态方程是直接描述系统输入和输出量之间的制约关系，是连续控制系统其他数学模型表达式的基础.( 错）4．global可以定义全局变量，全局变量的作用域是该MATLAB函数的整个工作区，其他的函数不能对它们进行存取和修改( 错）5．MATLAB中clf用于清除图形窗口上的旧图形（对）6.控制系统的数学模型有状态空间表达式，微分方程和积分方程( 错）7.仿真就是利用模型(物理模型或数学模型）代替实际系统进行实验和研究（对）8.离散相似法采样周期的选择应该满足香农定理(采用定理）（对）9.通常情况下，模拟仿真较数字仿真精度高( 错）10.机理模型法需要对系统的内部结构和特性完全的了解，但其精度较低( 错）11.绘制系统根轨迹的命令式是rlocus（对）12.仿真所遵循的基本原理是相似原理，即几何相似和数学相似（对）13.在MATLAB中，plot命令用于绘制三维图形( 错）14.绘制系统单位阶跃响应曲线的命令是step（对）15.系统仿真有三个基本的活动是模型建立，模型变换和模拟实验( 错）16.机理模型法就是对已知结构，参数的物理系统运用相应的物理定律或定理，经过合理的分析简化建立起来的各物理量间的关系（对）17.欧拉法的计算精度p=3( 错）18.绘制系统单位脉冲响应曲线的命令是implus（对）19.MATLAB的含义为矩阵实验室（对）20.margin(G)的含义是计算系统的相角裕度和幅值裕度（对）四、问答题:1.什么是仿真?它的主要优点是什么?它所遵循的基本原则是什么?[答案]:系统仿真是以相似原理，系统技术，信息技术及其应用领域有关的专业技术为基础，以计算机和各种专用物理效应设备为工具，利用系统模型对真实的或设想的系统进行动态研究的一门多学科的综合性技术.它是非常重要的设计自动控制系统或者评价系统性能和功能的一种技术手段.仿真的主要优点是:方便快捷，成本低廉，工作效率和计算精度都很高.它所遵循的基本原则是相似性原理.2.控制系统CAD可解决那些问题?[答案]:控制系统CAD可以解决以频域法为主要内容的经典控制理论和以时域法为主要内容的现代控制理论.此外，自适应控制，自校正控制以及最优控制等现代控制测略都可利用CAD 技术实现有效的分析与设计.3.控制系统建模的基本方法有哪些?他们的区别和特点是什么?[答案]:控制系统的建模方法大体有三种:机理模型法，统计模型法和混合模型法.机理模型法就是对已知结构，参数的物理系统运用相应的物理定律或定理，经过合理的分析简化建立起来的各物理量间的关系.该方法需要对系统的内部结构和特性完全的了解，精度高.统计模型法是采用归纳的方法，根据系统实测的数据，运用统计规律和系统辨识等理论建立的系统模型.该方法建立的数学模型受数据量不充分，数据精度不一致，数据处理方法的不完善，很难在精度上达到更高的要求.混合法是上述两种方法的结合.4.什么是离散系统?什么是离散事件系统?如何用数学的方法描述它们?[答案]:本课程所讲的”离散系统”指的是离散时间系统，即系统中状态变量的变化仅发生在一组离散时刻上的系统.它一般采用差分方程，离散状态方程和脉冲传递函数来描述.离散事件系统是系统中状态变量的改变是由离散时刻上所发生的事件所驱动的系统.这种系统的输入输出是随机发生的，一般采用概率模型来描述.5.动态系统仿真中常用的数值算法有哪几类，分别是什么?[答案]:主要有求解线性和非线性微分方程的数值积分法和计算线性时不变动态系统的离散相似法.其中，数值积分法主要有:欧拉(Euler)法，梯形法，龙格—库塔(Runge-Kutta)法和阿达姆斯(Adams)法；离散相似法主要有:置换法和相似变换法.6.为什么说模拟仿真较数字仿真精度低?其优点如何?.[答案]:由于受到电路元件精度的制约和容易受到外界的干扰，模拟仿真较数字仿真精度低，但模拟仿真具有如下优点:(1)描述连续的物理系统的动态过程比较自然和逼真；(2)仿真速度极快，失真小，结果可信度高；(3)能快速求解微分方程.模拟计算机运行时各运算器是并行工作的，模拟机的解题速度与原系统的复杂程度无关；(4)可以灵活设置仿真试验的时间标尺，既可以进行实时仿真，也可以进行非实时仿真；(5)易于和实物相连.7.采样控制系统数字仿真中连续部分离散化时的步长h如何选取?[答案]:(1)若仿真的任务仅要求计算系统输出y(t)而不要求计算系统内部状态变量，且连续部分的整体脉冲传递函数G(z)=Z[Gh(s)G0(s)]较易求出时，可选h=T(2)若连续部分整体脉冲传递函数G(z)=Z[Gh(s)G0(s)]不易求出；或仿真的任务要求计算系统输出y(t)和内部状态变量；或被控对象含有非线性环节时，可选h=T/N(N为正整数）.8.采样控制系统数字仿真有哪几种方法?[答案]:采样控制系统仿真通常有差分方程递推求解法，双重循环方法，应用MATLAB控制工具箱时域响应分析函数法和Simulink仿真法.9.计算机仿真有哪些优点?[答案]:(1）对尚处于论证或设计阶段的系统进行研究，唯一的方法就是仿真.(2）经济，安全，效率高.(3）研究系统非常方便灵活.10.评价优化方法的优劣的应该考虑哪些因素?[答案]:三方面因素:(1）收敛性:收敛性的好坏表示某种优化方法适用范围的大小，具体表示算法对于相当一类目标函数均能找到最优点.(2）收敛速度:为了求出同样精度的最优点，不同的优化方法所需要的迭代次数不同，迭代次数少的优化方法收敛速度较快.(3）每步迭代所需的计算量:每步迭代所需的计算量也是决定寻优速度的另一重要因素.。
控制系统数字仿真 四阶龙格库塔法
控制系统数字仿真1.实验目的1.掌握利用四阶龙格-库塔（Runge-Kutta）法进行控制系统数字仿真的方法。
2.学习分析高阶系统动态性能的方法。
3.学习系统参数改变对系统性能的影响。
二、实验内容已知系统结构如下图若输入为单位阶跃函数，计算当超调量分别为5%，25%，和50%时K的取值（用主导极点方法估算），并根据确定的K值在计算机上进行数字仿真。
三、实验过程1.计算K值二阶系统单位阶跃响应的超调量%100%=⨯1.当σ%=5%时解得 ζ=0.690设主导极点=ζa + a=0.69a+j0.72a代入D （s ）= 321025s s s K +++=0中， 32(0.690.72)10(0.690.72)25(0.690.72)0a j a a j a a j a K ++++++=解得K=31.3，a=-2.10即1,21.45 1.52s j =-±2. 当σ%=25%时解得 ζ=0.403设主导极点=ζa + a=0.403a+j0.915a代入D （s ）= 321025s s s K +++=0中， 32(0.4030.915)10(0.4030.915)25(0.4030.915)0a j a a j a a j a K ++++++=解得K=59.5，a=-2.75即1,21.112.53s j =-±3. 当σ%=50%时解得 ζ=0.215设主导极点=ζa + a=0.215a+j0.977a代入D （s ）= 321025s s s K +++=0中， 32(0.2150.977)10(0.2150.977)25(0.2150.977)0a j a a j a a j a K ++++++=解得K=103，a=-3.48即1,20.75 3.4s j =-±1. 计算调节时间和超调量 将不同K 值带入到程序中，利用四阶龙格-库塔法得到如下结果：1.K=31.3时， Ts=0.7550S, σ%=4.70% 2.K=59.5时， Ts=1.4100S ，σ%=23.28% 3.K=103时， Ts=1.9700S, σ%=45.49% 1. 用MATLAB 绘制2()(5)K G S S S =+的根轨迹图如下2. 绘制降阶系统跃响应曲线对原系统进行降阶处理，所得闭环传递函数为2()()1025C S K R S S S K=++, 利用四阶龙格-库塔法绘制阶跃响应曲线如下： -25-20-15-10-50510-15-10-551015Root LocusReal Axis I m a g i n a r y A x i s2.K=59.51.验证精确K值通过程序验证得到的精确K值分别为：K=31.76(σ%=5%);K=62.48(σ%=25%); K=113.82(σ%=50%)四、实验结论1.将系统传递函数化成时域形式，可以得到一组微分方程，利用四阶龙格-库塔法，就可以计算得到系统的响应。
四阶龙格-库塔解微分方程方法
计算机仿真作业2（院、系） 专业 班一、实验目的熟悉matlab 环境和基本操作。
二、实验内容熟悉matlab 环境及在自动控制中的应用。
面向微分方程的数字仿真1. 分别使用（1）四阶龙格-库塔解微分方程方法、（2）用matab 的ode45函数求解具有如下闭环传递函数的系统的阶跃响应。
43210()8364010s s s s s φ=++++解：（1）四阶龙格-库塔解微分方程方法r=1numo=[10];demo=[1 8 36 40 10];numh=0.002;demh=1;[num,den]=feedback(numo,demo,numh,demh);[A,B,C,D]=tf2ss(num,den);X=[zeros(length(A),1)];Y=0;t=0;Tf=5;h=0.02;n=Tf/h;for i=1:nK1=A*X+B*r;K2=A*(X+h*K1/2)+B*r;K3= A*(X+h*K2/2)+ B*r;K4=A*(X+h*K3)+B*r;X=X+h*(K1+2*K2+2*K3+K4)/6;Y=[Y,C*X+D*r];t=[t,t(i)+h];endplot(t,Y)(2) ode45函数解微分方程法：(2) ode45函数解微分方程法1、把连续函数转换成一阶微分方程组：用MATLAB编程>> num=[10];den=[1 8 36 40 10];[A B C D]=tf2ss(num,den)A =-8 -36 -40 -101 0 0 00 1 0 00 0 1 0B =1C =0 0 0 10D =0 所以模型为:function fun=fun(t,x)u=100;fun=[-8*x(1)-36*x(2)-40*x(3)-10*x(4)+u;x(1);x(2);x(3)]x0=[0,0,0,0],tspa=[0,5];x0 =0 0 0 0>> [t,y]=ode45('fun',tspa,x0)plot(t,y)[]x y u x x 10 0 0 000010 1 0 0 0 0 1 0 0 0 0 1 10- 40- 36- 8-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=。
四阶龙格库塔方法求解n自由度二阶微分方程
四阶龙格库塔方法求解n自由度二阶微分方程在数值计算中，龙格-库塔方法（Runge-Kutta method）是一种常用的求解常微分方程（ODE）的数值方法。
它是一种多步法（multiple-step method），通过使用不同的近似值来迭代计算ODE 的解。
其中，四阶龙格-库塔方法是最常用的龙格-库塔方法之一，也是一种高阶方法，可用于求解n自由度二阶微分方程。
下面将介绍四阶龙格-库塔方法的原理、步骤和应用。
首先，我们来回顾一下二阶微分方程的一般形式：y''(y)=y(y,y(y),y'(y))这里，y(y,y(y),y'(y))是已知的函数，y(y)是我们要求解的未知函数。
为了求解这个二阶微分方程，我们需要将其转化为一个一阶微分方程的求解问题。
令y(y)=y(y)y(y)=y'(y)这样，我们可以将原始的二阶微分方程转化为如下的一阶微分方程组：y'(y)=y(y)y'(y)=y(y,y(y),y(y))现在，我们可以利用四阶龙格-库塔方法来求解这个一阶微分方程组。
四阶龙格-库塔方法基于泰勒展开的思想，通过使用一系列的近似值来逼近方程的解。
该方法的一般形式为：y1=ℎy(yy,yy)y2=ℎy(yy+ℎ/2,yy+y1/2)y3=ℎy(yy+ℎ/2,yy+y2/2)y4=ℎy(yy+ℎ,yy+y3)yy+1=yy+1/6(y1+2y2+2y3+y4)yy+1=yy+ℎ其中，yy是当前时间步，yy+1是下一个时间步，yy是当前位置的近似解，yy+1是下一个位置的近似解，ℎ是时间步长，y1、y2、y3和y4是中间的近似值。
四阶龙格-库塔方法的步骤如下：Step 1:给定初始条件我们需要提供初始条件，包括初始位置y0和初始速度y0，以及时间步长ℎ、求解的时间区间[y0,yy]和离散时间步数y。
y0=y(y0)y0=y'(y0)Step 2:进行迭代计算从n=0开始，进行y步的迭代计算，其中y=(yy−y0)/ℎ。
四阶龙格——库塔法
四阶龙格——库塔法2013-2014（1）专业课程实践论文题目：四阶龙格—库塔法一、算法理论由定义可知，一种数值方法的精度与局部截断误差()po h有关，用一阶泰勒展开式近似函数得到欧拉方法，其局部截断误差为一阶泰勒余项2()o h，故是一阶方法，完全类似地若用p阶泰勒展开式2'''()11()()()......()()2!!pp p n n n n n h h y y x hy x y x y x O h p ++=+++++ 进行离散化，所得计算公式必为p 阶方法，式中'''''()(,),()(,)(,)(,)....x y x f x y y x f x y f x y f x y ==++由此，我们能够想到，通过提高泰勒展开式的阶数，可以得到高精度的数值方法，从理论上讲，只要微分方程的解()y x 充分光滑，泰勒展开方法可以构造任意的有限阶的计算公式，但事实上，具体构造这种公式往往相当困难，因为符合函数(,())f x y x 的高阶导数常常是很烦琐的，因此，泰勒展开方法一般不直接使用，但是我们可以间接使用泰勒展开方法，求得高精度的计算方法。
首先，我们对欧拉公式和改进欧拉公式的形式作进一步的分析。
如果将欧拉公式和改进的欧拉公式改写成如下的形式：欧拉公式{111(,)n n n n y y hK K f x y +==+改进的欧拉公式11211()22n n y y h K K +=++， 1(,)n n K f x y =，21(,)n n K f x h y hK =++。
这两组公式都是用函数(,)f x y 在某些点上的值的线性组合来计算1()n y x +的近似值1n y +，欧拉公式每前进一步，就计算一次(,)f x y 的值。
另一方面它是1()n y x +在n x 处的一阶泰勒展开式，因而是一阶方法。
改进的欧拉公式每前进一步，需要计算两次(,)f x y 的值。
最新matlab 四阶龙格-库塔法求微分方程
m a t l a b四阶龙格-库塔法求微分方程Matlab 实现四阶龙格-库塔发求解微分方程从理论上讲，只要函数在某区间上充分光滑，那么它可以展开为泰勒级数，因此在该区间上的函数值可用各阶导数值近似地表示出来，反之其各阶导数值也可用某些函数值的线性组合近似地表示出来。
龙格-库塔法就是将待求函数)(t y 展开为泰勒级数，并用方程函数),(y f t 近似其各阶导数，从而迭代得到)(t y 的数值解。
具体来说，四阶龙格-库塔迭代公式为)22(6143211k k k k h n n ++++=+y y ),(1n n t k y f =)2/,2/(12hk h t k n n ++=y f)2/,2/(23hk h t k n n ++=y f),(33hk h t k n n ++=y f实验内容：已知二阶系统21x x= ，u x x x 5.02.04.0212+--= ，0)0()0(21==x x ，u 为单位阶跃信号。
用四阶龙格-库塔法求数值解。
分析步长对结果的影响。
实验总结：实验报告要求简要的说明实验原理；简明扼要地总结实验内容；编制m 文件，并给出运行结果。
报告格式请按实验报告模板编写。
进入matlab ，Step1：choose way1 or way2way1）：可以选择直接加载M 文件（函数M 文件）。
way2）：点击new ——function ，先将shier （函数1文本文件）复制运行;点击new——function，再将RK（函数2文本文件）运行;点击new——function，再将finiRK（函数3文本文件）运行;Step2：回到command页面输入下面四句。
[t,k]=finiRK45([0;0],150);%迭代150次，步长=20/150[t1 k1]=ode45(@shier,[0 -10],[0 0]);%调用matlab自带四阶龙格-库塔，对比结果[t2 k2]=ode45(@shier,[0 10],[0 0]);plot(t,k(1,:),'-',t1,k1(:,1),'*',t2,k2(:,1),'^')%在图形上表示出来补充：改变步长影响数据的准确性。
四阶runge库塔法例题
四阶runge库塔法例题四阶Runge-Kutta法（也称为RK4法）是一种常用的数值积分方法，用于求解常微分方程的初值问题。
下面我将以一个例题来说明如何使用四阶Runge-Kutta法。
假设我们要求解如下的常微分方程初值问题：dy/dx = x^2 + y^2。
y(0) = 1。
首先，我们需要将该常微分方程转化为一阶形式。
令u = y，则dy/dx = du/dx。
将原方程代入可得：du/dx = x^2 + u^2。
u(0) = 1。
接下来，我们可以按照以下步骤使用四阶Runge-Kutta法进行数值求解：步骤1，确定积分步长h和积分区间。
假设我们要在区间[0, 1]上求解，可以选择一个适当的步长h，比如h = 0.1。
步骤2，初始化。
令x = 0，u = 1。
步骤3，根据四阶Runge-Kutta法的迭代公式进行迭代计算。
首先计算k1：k1 = h (x^2 + u^2)。
然后计算k2：k2 = h ((x + h/2)^2 + (u + k1/2)^2)。
接着计算k3：k3 = h ((x + h/2)^2 + (u + k2/2)^2)。
最后计算k4：k4 = h ((x + h)^2 + (u + k3)^2)。
步骤4，更新u和x的值。
u = u + (k1 + 2k2 + 2k3 + k4) / 6。
x = x + h.步骤5，重复步骤3和步骤4，直到x达到积分区间的上限。
通过以上步骤的迭代计算，我们可以得到在给定积分区间上的数值解。
需要注意的是，四阶Runge-Kutta法是一种较为精确的数值积分方法，但在某些特殊情况下可能会出现数值误差累积的问题。
因此，在实际应用中，需要根据具体问题的特点选择合适的数值积分方法。
希望以上回答能够满足你的要求。
如果还有其他问题，请随时提问。
实验：控制系统数字仿真之数值积分法
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
0 50 100 150 200 250 300 350 400 450 500
0
50
100
150
200
250
300
350
400
450
500
1.4
1.2
ts=199, Mp=18.2175, FAI=0.93537, tr=73, tp=106, ys=1.0003 ts=147, Mp=16.7351, FAI=0.94362, tr=74, tp=107, ys=1.0003 ts=147, Mp=16.7505, FAI=0.94355, tr=74, tp=107, ys=1.0003
ts=147.2, Mp=16.8911, FAI=0.94279, tr=73.9, tp=107, ys=1.0003
1
0.8
0.6
0.4
0.2
0
0
50
100
150
200
250
300
350
400
450
500
用梯形法得出系统响应曲线：
若采用欧拉法，误差为红色曲线围成的面积，而如果用梯形法，误差减少为 蓝色曲线围成的面积。同时，要求出蓝色曲线围成的面积，就要先出下一个点的 值。因此增加了计算量。 算法： 先用欧拉法求出下一个点的值， 用下一个点的值求这个点的斜率， 接着就能 求出梯形的面积。用新的面积（代表斜率）求出下一个点的值。 实验程序代码（与之前相同的部分没有复制）：
实验：控制系统数字 仿真之数值积分法
实验目的：
matlab用四阶龙格库塔函数求解微分方程组
一、介绍Matlab作为一种强大的科学计算软件，提供了众多函数和工具来解决微分方程组。
其中，四阶龙格库塔函数是一种常用的数值方法，用于求解常微分方程组。
本文将介绍如何使用Matlab中的四阶龙格库塔函数来求解微分方程组，并对该方法的原理和实现进行详细说明。
二、四阶龙格库塔方法四阶龙格库塔方法是一种常用的数值方法，用于求解常微分方程组。
它是一种显式的Runge-Kutta方法，通过逐步逼近微分方程的解，在每一步使用多个中间值来计算下一步的解。
该方法通过四个中间值来计算下一步的状态，并且具有较高的精度和稳定性。
三、在Matlab中使用四阶龙格库塔方法求解微分方程组在Matlab中，可以使用ode45函数来调用四阶龙格库塔方法来解决微分方程组的问题。
ode45函数是Matlab提供的用于求解常微分方程组的函数，可以通过指定微分方程组以及初值条件来调用四阶龙格库塔方法来进行求解。
1. 定义微分方程组我们需要定义要求解的微分方程组。
可以使用Matlab中的匿名函数来定义微分方程组，例如：```matlabf = (t, y) [y(2); -sin(y(1))];```其中，f是一个匿名函数，用于表示微分方程组。
在这个例子中，微分方程组是y' = y2, y2' = -sin(y1)。
2. 指定初值条件和求解区间接下来，我们需要指定微分方程组的初值条件和求解区间。
初值条件可以通过指定一个初始时刻的状态向量来完成，例如：```matlabtspan = [0, 10];y0 = [0, 1];```其中，tspan表示求解区间，y0表示初值条件。
3. 调用ode45函数进行求解我们可以通过调用ode45函数来求解微分方程组的数值解。
具体的调用方式如下：```matlab[t, y] = ode45(f, tspan, y0);```其中，t和y分别表示求解的时间点和对应的状态值。
四、示例下面我们通过一个具体的例子来演示如何使用Matlab中的四阶龙格库塔方法来求解微分方程组。
四阶龙格库塔法
龙格库塔方法二：
当系统的稳定在一定范围内时，龙格库塔法可采用以下方法：
其中
要给定一个特定的方法，必须提供级数s以及系数aijbi和ci.龙格库塔表如下：
0
龙格库塔法是自治的，如果
对第一种方法，龙格库塔表如下：
0
1/2
经典四阶龙格库塔法
系统方程和表述如下：
则系统的输出按如下求解：
其中：
这样，下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均：
k1是时间段开始时的斜率；
k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn + h/2的值；
k3也是中点的斜率，但是这次采用斜率k2决定y值；
1/2
1/2
0
1/2
1
0
0
1
1/6
1/3
1/3
1/6
龙格库塔方法三：
对刚性系统，龙格库塔方法具有以下形式：
其中：
就得到四阶龙格库塔法
1.00000000000000 0.96078943915232 0.85214378896621 0.69767632607103 0.52729242404305 0.36787944117144 0.23692775868212
2-3-3 梯形法（RK2） 用梯形面积代替曲边梯形的面积。
2-3 数值积分法的稳定性 三、用测试方程考察欧拉法的稳定性 1. xn 1 xn hxn ) (1 h ) xn (1 h ) n 1 x0
当|1+h |<=1时，计算稳定。 2、求其稳定边界 设h =x+jy, 由|1+h |<=1 得:|1+x+jy|<=1 稳定边界为：(1+x)2+y2=1 当为负实数时，步长的 稳定区间为：0<h<|2/ |
我们所讨论的系统主要以线性定常连续系统为 主。 线性系统满足叠加性和齐次性。
1. 微分方程建立的一般步骤 采用解析法来建立系统或元部件的微分方程所遵循的一般步 骤是： （1）确定系统或元部件的输入、输出变量。 （2）根据物理和化学定律（比如：牛顿运动定律、能量守 恒定律、克希霍夫定律等）列出系统或元部件的原始方程式， 按照工作条件忽略一些次要因素。 （3）找出原始方程式中间变量与其它因素的关系式。 （4）消去原始方程式的中间变量，得到一个关于输入、输 出的微分方程式。 （5）进行标准化处理，将输出各项放在等号左端，输入各 项放在等号右端，并且按照微分方程的阶次降幂排列，同时 将各系数化为具有一定物理意义的形式。
3、状态空间表达式
AX BU X Y CX DU
一阶微分方程组
2.1.2 描述控制系统常用的数学模型
四阶龙格库塔法原理
四阶龙格库塔法原理
四阶龙格库塔法是一种常见的数值计算方法，用于求解常微分方程的初值问题。
该方法通过逐步逼近准确解来得到数值解。
在四阶龙格库塔法中，我们将求解区间 [a, b] 平均分成 n 个子
区间，每个子区间的长度为 h = (b - a) / n，其中 a 是初始点，b 是终点。
首先，我们需要给出初始条件 y(a)，即方程在 a 点的值。
然后，我们利用以下公式依次计算 y(a + kh) 的近似值，其中 k 为当
前步数（从 0 开始）：
\[
\begin{align*}
k_1 &= hf(a, y(a)) \\
k_2 &= hf(a + \frac{h}{2}, y(a) + \frac{k_1}{2}) \\
k_3 &= hf(a + \frac{h}{2}, y(a) + \frac{k_2}{2}) \\
k_4 &= hf(a + h, y(a) + k_3) \\
y(a + kh) &= y(a) + \frac{k_1 + 2k_2 + 2k_3 + k_4}{6}
\end{align*}
\]
其中，f(a, y(a)) 为方程的右端函数，描述了方程随时间的变化
规律。
在每一步中，我们利用当前步数的近似值 k1 到 k4 来计算下一步的近似值。
最后，通过不断迭代计算，直到达到指定的终点b，我们可以得到方程的数值解。
四阶龙格库塔法的主要原理是通过将步长 h 下的误差控制在O(h^5) 的级别，从而在保证计算效率的同时，提高数值解的精度。
该方法具有较好的稳定性和精确性，常被应用于各种科学计算和工程问题的数值求解中。
定步长四阶龙格-库塔(Runge-Kutta)法
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''' 模块名：定步长四阶龙格_库塔_Runge_Kutta_法.bas' 功能：解一阶常微分方程组的初值问题。
子过程Runge_Kutta由给定步长h和初始点上的值用四阶龙格-库塔法' 求解给定的初值问题，调用它一次向前积分一步，一般用于求解某一小区间中某几点的函数值；' 子过程Runge_KuttaDumb是连续调用于过程Runge_Kutta()计算出某一区间上的值，此时该区间长度可较' 大，但积分步长都较小。
一般来说，在精度要求不高时，可采用该方法''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''Option Explicit''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''' 模块名：定步长四阶龙格_库塔_Runge_Kutta_法.bas' 函数名：Runge_KuttaDumb(n，x1，x2，nstep，vstart()，xx()，y())' 功能：解一阶常微分方程组的初值问题。
matlab经典的4级4阶runge kutta法 -回复
matlab经典的4级4阶runge kutta法-回复什么是Matlab经典的4级4阶Runge-Kutta法（RK4）？Matlab经典的4级4阶Runge-Kutta法是一种数值解法，用于求解常微分方程（ODEs）。
该方法是由德国数学家卡尔·雷特克（Carl Runge）和瓦尔特·库塔（Martin Kutta）在1900年至1901年期间独立发现和发展的。
它是一种四阶精度的方法，具有较高的稳定性和所需步骤数较少的优点，因此被广泛应用于科学和工程领域的数值模拟和计算中。
Runge-Kutta法的基本思想是通过逐步逼近来计算给定的常微分方程。
根据初值条件，我们可以设置一个初始点并迭代地计算出下一个点的近似值。
使用RK4方法，我们可以通过四个步骤来计算下一个点的近似值。
那么，RK4的四个步骤是什么呢？首先，我们设定初值条件。
对于一个一阶常微分方程y' = f(x, y)，我们需要给出初始点(x0, y0)。
然后，我们选择一个适当的步长h，以确定我们在每个步骤中应该前进的距离。
其次，我们计算k1，这是在初始点上斜率的近似值。
我们使用初始点的坐标(x0, y0)和方程f(x, y)来计算k1 = f(x0, y0)。
接下来，我们计算k2，这是在前进到下一个点之前我们在中间点上斜率的近似值。
中间点的坐标为(x0 + h/2, y0 + h*k1/2)。
我们使用方程f(x, y)和中间点的坐标来计算k2 = f(x0 + h/2, y0 + h*k1/2)。
然后，我们计算k3，这是在前进到下一个点之前我们在另一个中间点上斜率的近似值。
另一个中间点的坐标为(x0 + h/2, y0 + h*k2/2)。
我们使用方程f(x, y)和另一个中间点的坐标来计算k3 = f(x0 + h/2, y0 + h*k2/2)。
最后，我们计算k4，这是在下一个点上的斜率的近似值。
下一个点的坐标为(x0 + h, y0 + h*k3)。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
控制系统数字仿真
1.实验目的
1.掌握利用四阶龙格-库塔（Runge-Kutta）法进行控制系统数字仿真的方
法。
2.学习分析高阶系统动态性能的方法。
3.学习系统参数改变对系统性能的影响。
二、实验内容
已知系统结构如下图
若输入为单位阶跃函数，计算当超调量分别为5%，25%，和50%时K的取值（用主导极点方法估算），并根据确定的K值在计算机上进行数字仿真。
三、实验过程
1.计算K值
二阶系统单位阶跃响应的超调量
%100%
=⨯
1.当σ%=5%时
解得 ζ=0.690
设主导极点
=ζa + a=0.69a+j0.72a
代入D （s ）= 32
1025s s s K +++=0中， 32(0.690.72)10(0.690.72)25(0.690.72)0
a j a a j a a j a K ++++++=解得K=31.3，a=-2.10
即1,2
1.45 1.52s j =-±
2. 当σ%=25%时
解得 ζ=0.403
设主导极点
=ζa + a=0.403a+j0.915a
代入D （s ）= 321025s s s K +++=0中， 32(0.4030.915)10(0.4030.915)25(0.4030.915)0
a j a a j a a j a K ++++++=解得K=59.5，a=-2.75
即1,2
1.11
2.53s j =-±
3. 当σ%=50%时
解得 ζ=0.215
设主导极点
=ζa + a=0.215a+j0.977a
代入D （s ）= 321025s s s K +++=0中， 32(0.2150.977)10(0.2150.977)25(0.2150.977)0
a j a a j a a j a K ++++++=解得K=103，a=-3.48
即1,2
0.75 3.4s j =-±
1. 计算调节时间和超调量 将不同K 值带入到程序中，利用四阶龙格-库塔法得到如下结果：
1.
K=31.3时， Ts=0.7550S, σ%=4.70% 2.
K=59.5时， Ts=1.4100S ，σ%=23.28% 3.
K=103时， Ts=1.9700S, σ%=45.49% 1. 用MATLAB 绘制2()(5)
K G S S S =+的根轨迹图如下
2. 绘制降阶系统跃响应曲线
对原系统进行降阶处理，所得闭环传递函数为
2()()1025C S K R S S S K
=++, 利用四阶龙格-库塔法绘制阶跃响应曲线如下： -25-20-15-10
-50510-15-10
-5
5
10
15
Root Locus
Real Axis I m a g i n a r y A x i s
2.K=59.5
1.验证精确K值
通过程序验证得到的精确K值分别为：K=31.76(σ%=5%);
K=62.48(σ%=25%); K=113.82(σ%=50%)
四、实验结论
1.将系统传递函数化成时域形式，可以得到一组微分方程，利用四阶龙格-库塔法，就可以计算得到系统的响应。
当然，这是一种近似解。
2.利用主导极点法，可以将高阶系统进行降阶，用二阶系统近似来分析。
3.开环系统的参数对闭环系统动态性能造成影响：当开环比例系数适当，系统动态性能较好的情况下，用主导极点的方法，不至于造成较大的误差；当开环比例系数较大，系统动态性能较差时，采取同样的方法，产生了较大的误差。
程序清单
A=[0 1 0;0 0 1;-k -25 -10];
b=[0 0 1]';
c=[k 0 0];
X=zeros(3,1);
t=0:0.01:10;
n=length(t);
h=0.01;
for i=1:n
K1=A*X+b;
K2=A*(X+(h/2)*K1)+b;
K3=A*(X+(h/2)*K2)+b;
K4=A*(X+h*K3)+b;
X=X+(h/6)*(K1+2*K2+2*K3+K4);
y(i)=c*X;
end
plot(y);
s=1001;while y(s)>0.95&y(s)<1.05;s=s-1;end;
t=(s-1)*0.005；max(y)-1。