自动控制原理实验四-控制系统数字仿真四阶龙格库塔法

合集下载

四阶龙格库塔实验报告

三、四阶Runge-Kutta 法求解常微分方程一、龙格库塔法的思想根据第九章的知识可知道，Euler 方法的局部截断误差是2()O h ，而当用Euler 方法估计出1,()(1)n n n n y y hf x y +=+ 再用梯形公式111[(,)(,)](2)2n n n n n n h y y f x y f x y +++=++进行校正，即采用改进Euler 方法得出数值解的截断误差为3()O h 。

由Lagrange 微分中值定理'11()()()()()(,())(3)n n n n n y x y x y x x y x hf y ξξξ++=+-=+ 记*(,())k hf y ξξ=，得到*1()()(4)n n y x y x k +=+这样只要给出一种计算*k 的算法，就能得到相应的计算公式。

用这种观点的来分析Euler 方法和改进Euler 方法，Euler 方法的迭代公式可改写为111(,)n n n n y y k k hf x y +=+=改进Euler 方法的预报-校正公式可改写为 1121211()2(,),(,)n n n n n n y y k k k hf x y k hf x h y k +=++==++ Euler 方法实际上是用一个点处的值1k 近似*k ，而改进Euler 方法是用两个点处的值1k ，和2k ，做算术平均值近似*k 自然改进Euler 方法要优于Euler 方法。

因此，可以想到假如在1[,]n n x x +内多预报几个点值i k ，并用他们的加权平均值作为*k 的近似值，则有可能构造出具有更高精度的计算公式，这就是Runge-Kutta 法的基本思想。

二、四阶龙格库塔法由Runge-Kutta 的基本思想，构造四阶Runge-Kutta 法是利用1234,,k k k k 和的加权平均值来近似*k ，因此令1112233441211132211243312213(,)(,)(5)(,)(,)n n n n n n n n n n y y w K w K w K w K K hf x y K hf x h y K K hf x h y K K K hf x h y K K K αβαβγαβγη+=++++⎧⎪=⎪⎪=++⎨⎪=+++⎪⎪=++++⎩使得511()()n n y x y O h ++-=即其总体截断误差为4()O h 。

四阶龙格-库塔法求解常微分方程的初值问题-matlab通用程序

参考教材《数值分析》李乃成.梅立泉clearclcformat longm=input('请输入常微分方程的阶数m=');a=input('请输入x下限a=');b=input('请输入x上限b=');h=input('请输入步长h=');ym=input('令y(1,1)=y,y(2,1)=y’,y(3,1)=y’’...请输入ym=','s'); %输入的时候必须按照这个形式输入y1=y(1,1);if m==1 %一阶初值问题单独求解mm=(b-a)/h;y(1,1)=input('请输入在初值点的函数值f(a)=');x=a;y11(1)=y(1,1);for k1=2:(mm+1)y1=y(1,1);K(1,1)=h*(eval(ym)); %计算K1x=x+h/2;y(1,1)=y1+K(1,1)/2;y1=y(1,1);K(1,2)=h*(eval(ym)); %计算K2x=x;y(1,1)=y1+K(1,2)/2-K(1,1)/2;y1=y(1,1);K(1,3)=h*(eval(ym)); %计算K3x=x+h/2;y(1,1)=y1+K(1,3)-K(1,2)/2;y1=y(1,1);K(1,4)=h*(eval(ym)); %计算K4y11(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4))/6; y(1,1)=y11(k1);x=a+(k1-1)*h;endy11else %高阶初值问题mm=(b-a)/h; %一共要求解mm个数据点for k2=1:m %读取初值条件fprintf('请输入%d阶导数的初值f(%d)(a)=\n',(k2-1),(k2-1));y(k2,1)=input('=');endfor k2=1:my22(1,k2)=y(k2,1); %先把初值保存在矩阵y22(m,n)中,m表示第几个所求点,n表示第n阶初值endx=a;for k4=2:(mm+1) %求解mm个数据点的循环for k=1:(m-1) %计算K1,包括每一阶的K1 K(k,1)=h*y(k+1,1); %y(k+1,1)中k+1表示第k+1阶,1表示第一个点;K(k,1)中k表示阶数,1表示K1endK(m,1)=h*(eval(ym));x=x+h/2; %求解K1之前,先重新对x和y赋值for k3=1:my(k3,1)=y(k3,1)+K(k3,1)/2;endfor k=1:(m-1) %计算K2K(k,2)=h*y(k+1,1);endK(m,2)=h*(eval(ym));x=x;for k3=1:my(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;endfor k=1:(m-1) %计算K3K(k,3)=h*y(k+1,1);endK(m,3)=h*(eval(ym));x=x+h/2;for k3=1:my(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2; %这里容易出错endfor k=1:(m-1) %计算K4K(k,4)=h*y(k+1,1);endK(m,4)=h*(eval(ym));for k5=1:my22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4))/6; %这里,除了要求出下一个点的数值,还要求出相应的导数值endfor k6=1:m %除了对y(1,1)重新赋值外,还要对y(2,1)等重新赋值y(k6,1)=y22(k4,k6);endx=a+(k4-1)*h;endy22(:,1) end。

杨源实验四连续系统的数字仿真

实验四连续系统的数字仿真
一、实验目的
1、掌握四阶龙格-库塔法系统仿真的方法；
2、掌握面向结构图仿真的方法与步骤；
3、熟悉仿真程序的编写方法；
二、实验内容
1、SISO系统的传递框图如图所示，试编写程序用四阶龙格-库塔法对系统进行仿真，求系统输出y的动态响应
解：
输入命令：simulink
在新打开的窗口中选择File- ->New- ->Model
在untitled窗口中绘制图形如下：
杨源
单击仿真按钮Start simulation，然后双击示波器（Scope）显示如下结果
2、某一系统由四个典型环节组成，如图所示，求输出量y的动态响应。

解：
输入命令：simulink
在新打开的窗口中选择File- ->New- ->Model
在untitled窗口中绘制图形如下：
单击仿真按钮Start simulation，然后双击示波器（Scope）显示如下结果：。

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.算法原理对于一阶常微分方程组的初值问题⎪⎪⎪⎩⎪⎪⎪⎨⎧=⋯⋯==⋯⋯=⋯⋯⋯⋯=⋯⋯=0020********'212'2211'1)(,,)(,)())(,),(),(,()())(,),(),(,()())(,),(),(,()(n n n n n n n y x y y x y y x y x y x y x y x f x y x y x y x y x f x y x y x y x y x f x y ，其中b x a ≤≤。

若记Tn Tn Tn y x f y x f y x f y x f y y y y x y x y x y y x y )),(,),,(),,((),(),,,())(),(),(()(2102010021⋯⋯=⋯⋯=⋯⋯=，，则可将微分方程组写成向量形式⎩⎨⎧=≤≤=0')()),(,()(y a y b x a x y x f x y微分方程组初值问题在形式上和单个微分方程处置问题完全相同，只是数量函数在此变成了向量函数。

因此建立的单个一阶微分方程初值问题的数值解法，可以完全平移到求解一阶微分方程组的初值问题中，只不过是将单个方程中的函数转向向量函数即可。

标准4阶R-K 法的向量形式如下：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++=++=++==++++=+),()21,2()21,2(),()22(61342312143211K y h x hf K K y h x hf K K y h x hf K y x hf K K K K K y y n n n n n n n n n n 其分量形式为n j K y K y K y h x hf K K y K y K y h x hf K K y K y K y h x hf K y y y x hf K K K K K y y n ni i i i j j n nii i i j j n nii i i j j ni i i i j j j j j j i j i j ,,2,1).,,,;(),2,2,2;2(),2,2,2;2(),,,,;(),22(6132321314222212131212111221143211,1,⋯⋯=⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧+⋯⋯+++=+⋯⋯+++=+⋯⋯+++=⋯⋯=++++=++，，2.程序框图3.源代码%该函数为四阶龙格-库塔法function [x,y]=method(df,xspan,y0,h)%df为常微分方程，xspan为取值区间，y0为初值向量，h为步长x=xspan(1):h:xspan(2);m=length(y0);n=length(x);y=zeros(m,n);y(:,1)=y0(:);for i=1:n-1k1=feval(df,x(i),y(:,i));k2=feval(df,x(i)+h/2,y(:,i)+h*k1/2);k3=feval(df,x(i)+h/2,y(:,i)+h*k2/2);k4=feval(df,x(i)+h,y(:,i)+h*k3);y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;end%习题9.2clear;xspan=[0,1];%取值区间h=0.05;%步长y0=[-1,3,2];%初值df=@(x,y)[y(2);y(3);y(3)+y(2)-y(1)+2*x-3];[xt,y]=method(df,xspan,y0,h)syms t;yp=t*exp(t)+2*t-1;%微分方程的解析解yp1=xt.*exp(xt)+2*xt-1%计算区间内取值点上的精确解[xt',y(1,:)',yp1']%y(1,:)为数值解，yp1为精确解ezplot(yp,[0,1]);%画出解析解的图像hold on;plot(xt,y(1,:),'r');%画出数值解的图像4.计算结果。

控制系统数字仿真

底盘控制
对汽车的悬挂、转向、制动等系统进行数字仿真，验证底盘控制算法的正确性和可行性，提高汽车的操控稳定性和行驶安全性。
自动驾驶控制
通过数字仿真技术，模拟自动驾驶系统的行为和性能，评估自动驾驶控制算法的优劣和适用性，推动自动驾驶技术的发展和应用。
04
控制系统数字仿真挑战与解决方案
实时性挑战与解决方案
电机控制
对电机的启动、调速、制动等过程进行数字仿真，验证电机控制算法的正确性和可行性，提高电机的稳定性和可靠性。
智能控制
通过数字仿真技术，模拟智能控制系统的行为和性能，评估智能控制算法的优劣和适用性。
机器人控制
1 2 3
运动控制
对机器人的关节和末端执行器进行数字仿真，模拟机器人的运动轨迹和姿态，验证运动控制算法的正确性和可行性。
实时性挑战
在控制系统数字仿真中，实时性是一个关键的挑战。由于仿真过程中需要不断进行计算和控制，如果仿真时间过长，会导致控制延迟，影响系统的实时响应。
解决方案
为了解决实时性挑战，可以采用高效的算法和计算方法，如并行计算、分布式计算等，以提高仿真速度。同时，可以通过优化仿真模型和减少不必要的计算来降低仿真时间。
特点
数字仿真具有高效、灵活、可重复性等优点，可以模拟各种实际工况和参数条件，为控制系统设计、优化和故障诊断提供有力支持。
数字仿真的重要性
验证设计
通过数字仿真可以对控制系统设计进行验证，确保系统性能符合预期要求。
优化设计
数字仿真可以帮助发现系统设计中的潜在问题，优化系统参数和性能。
故障诊断
THANபைடு நூலகம்S
感谢观看
发展趋势
目前，数字仿真正朝着实时仿真、高精度建模、智能化分析等方向发展，为控制系统的研究和应用提供更强大的支持。

控制系统数字仿真四阶龙格库塔法

控制系统数字仿真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.将系统传递函数化成时域形式，可以得到一组微分方程，利用四阶龙格-库塔法，就可以计算得到系统的响应。

4阶Runge-Kutta法求解一阶常微分方程

4阶Runge-Kutta法求解一阶常微分方程以下是为大家整理的4阶Runge-Kutta法求解一阶常微分方程的相关范文，本文关键词为Runge-Kutta,求解,一阶,微分方程,mATLAb语，您可以从右上方搜索框检索更多相关文章，如果您觉得有用，请继续关注我们并推荐给您的好友，您可以在教育文库中查看更多范文。

《mATLAb语言及应用》大作业姓学学班名号院级：：：：题目编号：20XX年10月134阶Runge-Kutta法求解一阶常微分方程。

一、Runge-Kutta法的数学理论龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。

由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。

该算法是构建在数学支持的基础之上的。

龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。

如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。

一阶常微分方程可以写作：y'=f(x,y),使用差分概念。

(Yn+1-Yn)/h=f(xn,Yn)推出（近似等于，极限为Yn'）Yn+1=Yn+h*f(xn,Yn)另外根据微分中值定理，存在0这里K＝f(xn+th,Y(xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。

利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为o(h^5)的四阶龙格库塔公式：K1＝f(xn,Yn);K2=f(xn+h/2,Yn+(h/2)*K1);K3=f(xn+h/2,Yn+(h/2)*K2);K4=f(xn+h,Yn+h*K3);Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6)二、Runge-Kutta的算法和流程图在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为0(h5),被广泛应用于解微分方程的初值问题。

四阶龙格-库塔法求解常微分方程的初值问题

实验 4 四阶龙格-库塔法求解常微分方程的初值问题
1. 算法原理龙格—库塔法是一种求其准确解 y( x) 在一系列点 xi 处 y ( xi ) 的近似值 yi 的方法， yi 称为数值解。经典的四阶龙格—库塔公式为：
K1 hf ( xi , yi ) K 2 hf ( xi h , yi 1 K1 ) 2 2 h 1 K 3 hf ( xi , yi K 2 ) 2 2 K 4 hf ( xi h, yi K 3 ) yi 1 yi 1 ( K1 2 K 2 2 K 3 K 4 ) 6
2. 程序框图
开始
输入函数 f ，初始点 x0，初始向量 y0，步长 h，矩阵的等分数 N
N
xi = x0 ih
计算 K1，K2，K3，
i=i+1 否
K4，yi+1
i=N
是
输出一系列点 xi 和对应的数值解 yi
结束
3. 程序使用说明本程序使用 MATLAB 利用四阶龙格—库塔法来求解常微分方程的初值问题。源程序文件 “RK.m”为龙格—库塔法的源程序，x为一系列点 xi = x0 ih 组成的向量，y为数值解 yi 组成的矩阵，f为1阶微分方程的函数，x0为初始点，y0为初始向量（列向量），h为步长，N为矩阵的等分数。输入1阶微分方程的函数f，初始点x0，初始向量y0，步长h，矩阵的等分数N 后，在命令窗口输入[x,y]=RK(f,x0,y0,h,N)，回车后即可算出一系列点 xi 和对应的数值解 yi 。源程序文件 “RK2.m”是计算实习 9.2 算例的程序，直接运行后，在命令窗口输入 t1，t2，回车后即可得到算例结果，其中 t1 表示 y (1) 的近似值，t2 表示 y (1) 的精确值。 4. 算例计算结果此算例为课本 304 页计算实习 9.2. 将算例的高阶常微分方程转化为如下一阶微分方程组。

四阶龙格-库塔解微分方程方法

计算机仿真作业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-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=。

四阶龙格——库塔法

四阶龙格——库塔法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 的值。

控制系统仿真实验报告

控制系统仿真实验报告班级：测控1402班姓名：王玮学号：072018年01月实验一经典的连续系统仿真建模方法一实验目的:1 了解和掌握利用仿真技术对控制系统进行分析的原理和步骤。

2 掌握机理分析建模方法。

3 深入理解阶常微分方程组数值积分解法的原理和程序结构，学习用Matlab编写数值积分法仿真程序。

4 掌握和理解四阶Runge-Kutta法，加深理解仿真步长与算法稳定性的关系。

二实验内容:1. 编写四阶 Runge_Kutta 公式的计算程序，对非线性模型（3）式进行仿真。

（1）将阀位u 增大10％和减小10％，观察响应曲线的形状;（2）研究仿真步长对稳定性的影响，仿真步长取多大时RK4 算法变得不稳定（3）利用 MATLAB 中的ode45()函数进行求解，比较与（1）中的仿真结果有何区别。

2. 编写四阶 Runge_Kutta 公式的计算程序，对线性状态方程（18）式进行仿真（1）将阀位增大10％和减小10％，观察响应曲线的形状;（2）研究仿真步长对稳定性的影响，仿真步长取多大时RK4 算法变得不稳定（4）阀位增大10％和减小10％，利用MATLAB 中的ode45()函数进行求解阶跃响应，比较与（1）中的仿真结果有何区别。

三程序代码:龙格库塔:%RK4文件clccloseH=[,]';u=; h=1;TT=[];XX=[];for i=1:h:200k1=f(H,u);k2=f(H+h*k1/2,u);k3=f(H+h*k2/2,u);k4=f(H+h*k3,u);H=H+h*(k1+2*k2+2*k3+k4)/6;TT=[TT i];XX=[XX H];end;hold onplot(TT,XX(1,:),'--',TT,XX(2,:)); xlabel('time')ylabel('H')gtext('H1')gtext('H2')hold on水箱模型:function dH=f(H,u)k=;u=;Qd=;A=2;a1=;a2=;dH=zeros(2,1);dH(1)=1/A*(k*u+Qd-a1*sqrt(H(1)));dH(2)=1/A*(a1*sqrt(H(1))-a2*sqrt(H(2)));2编写四阶 Runge_Kutta 公式的计算程序，对线性状态方程（18）式进行仿真：1 阀值u对仿真结果的影响U=;h=1; U=;h=1;U=;h=1;2 步长h对仿真结果的影响:U=;h=5; U=;h=20;U=;h=39 U=;h=50由以上结果知,仿真步长越大,仿真结果越不稳定。

实验：控制系统数字仿真之数值积分法

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
用梯形法得出系统响应曲线：
若采用欧拉法，误差为红色曲线围成的面积，而如果用梯形法，误差减少为蓝色曲线围成的面积。同时，要求出蓝色曲线围成的面积，就要先出下一个点的值。因此增加了计算量。算法：先用欧拉法求出下一个点的值，用下一个点的值求这个点的斜率，接着就能求出梯形的面积。用新的面积（代表斜率）求出下一个点的值。实验程序代码（与之前相同的部分没有复制）：
实验：控制系统数字仿真之数值积分法
实验目的：

四阶龙格库塔法

k4是时间段终点的斜率，其y值用k值：
龙格库塔方法二：
当系统的稳定在一定范围内时，龙格库塔法可采用以下方法：
其中
要给定一个特定的方法，必须提供级数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
龙格库塔方法三：
对刚性系统，龙格库塔方法具有以下形式：
其中：

龙格库塔法-原理及程序实现

龙格库塔法一、基本原理：“龙格-库塔(Runge-Kutta)方法”是一种在工程上应用广泛的“高精度单步算法”。

由于此算法精度高，采取措施对误差进行抑制，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法，即：yi＋１＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6K1=f(xi,yi)K2=f(xi+h/2,yi+h*K1/2)K3=f(xi+h/2,yi+h*K2/2)K4=f(xi+h,yi+h*K3)通常所说的龙格-库塔法就是指四阶——龙格库塔法，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式。

（1）计算公式（1）的局部截断误差是。

龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算在一些点上的值，如四阶龙格-库塔法每计算一步需要计算四次的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”。

二、小程序#include<stdio.h>#include<math.h>#define f(x,y) (-1*(x)*(y)*(y))void main(void){double a,b,x0,y0,k1,k2,k3,k4,h;int n,i;printf("input a,b,x0,y0,n:");scanf("%lf%lf%lf%lf%d",&a,&b,&x0,&y0,&n); printf("x0\ty0\tk1\tk2\tk3\tk4\n");for(h=(b-a)/n,i=0;i!=n;i++){k1=f(x0,y0);k2=f(x0+h/2,y0+k1*h/2);k3=f(x0+h/2,y0+k2*h/2);k4=f(x0+h,y0+h*k3);printf("%lf\t%lf\t",x0,y0);printf("%lf\t%lf\t",k1,k2);printf("%lf\t%lf\n",k3,k4);y0+=h*(k1+2*k2+2*k3+k4)/6;x0+=h;}printf("xn=%lf\tyn=%lf\n",x0,y0);}运行结果：input a,b,x0,y0,n:0 5 0 2 20x0y0k1k2k3 k40.000000 2.000000-0.000000-0.500000-0.469238-0.8861310.250000 1.882308-0.885771-1.176945-1.129082-1.2800600.500000 1.599896-1.279834-1.295851-1.292250-1.2227280.750000 1.279948-1.228700-1.110102-1.139515-0.9901621.000000 1.000027-1.000054-0.861368-0.895837-0.7528521.2500000.780556-0.761584-0.645858-0.673410-0.5621891.5000000.615459-0.568185-0.481668-0.500993-0.4205371.7500000.492374-0.424257-0.361915-0.374868-0.3178552.0000000.400054-0.320087-0.275466-0.284067-0.2435982.2500000.329940-0.244935-0.212786-0.218538-0.1894822.5000000.275895-0.190295-0.166841-0.170744-0.1495632.7500000.233602-0.150068-0.132704-0.135399-0.1197033.0000000.200020-0.120024-0.106973-0.108868-0.0970483.2500000.172989-0.097256-0.087300-0.088657-0.0796183.5000000.150956-0.079757-0.072054-0.073042-0.0660303.7500000.132790-0.066124-0.060087-0.060818-0.0553054.0000000.117655-0.055371-0.050580-0.051129-0.0467434.2500000.104924-0.046789-0.042945-0.043363-0.0398334.5000000.094123-0.039866-0.036750-0.037072-0.0342024.7500000.084885-0.034226-0.031675-0.031926-0.029571xn=5.000000yn=0.076927。

四阶龙格—库塔法在电机数字仿真中的应用问题

四阶龙格—库塔法在电机数字仿真中的应用问题
孟传富;张明玉
【期刊名称】《山东工业大学学报》
【年(卷),期】1989(019)003
【摘要】探讨应用四阶龙格-库塔法进行电机数字仿真时常遇到的有关问题,包括对电机方程组病态性的判断和仿真时最佳步长的确定等,实践证明,应用本文提出的判据和方法,可使电机的数字仿真得到快速、准确的效果。

【总页数】4页(P75-78)
【作者】孟传富;张明玉
【作者单位】不详;不详
【正文语种】中文
【中图分类】TM343
【相关文献】
1.四阶龙格-库塔法在捷联惯导系统姿态解算中的应用 [J], 张春慧;吴简彤;何昆鹏;周雪梅
2.四阶龙格—库塔法的原理及其应用 [J], 冯建强;孙诗一;
3.四阶龙格-库塔法在捷联惯导系统姿态解算中的应用 [J], 张春慧;吴简彤;何昆鹏;郭新伟
4.四阶龙格—库塔法在电机数字仿真中的适用范围与最佳步长的确定 [J], 邬学义
5.四阶龙格_库塔法在火控解算中的应用 [J], 李丹
因版权原因，仅展示原文概要，查看原文内容请购买。

北航自控实验报告

成绩北京航空航天大学自动控制原理实验报告学院仪器科学与光电工程学院专业方向测控班级学号学生姓名指导教师自动控制与测试教学实验中心实验一一、二阶系统的电子模拟及时域响应的动态测试实验时间 2013.10.22实验编号1 同组同学无一、实验目的1. 了解一、二阶系统阶跃响应及其性能指标与系统参数之间的关系。

2. 学习在电子模拟机上建立典型环节系统模型的方法。

3. 学习阶跃响应的测试方法。

二、实验内容1. 建立一阶系统的电子模型,观测并记录在不同时间常数t时的跃响应曲线,并测定其过渡过程时间ts。

2. 建立二阶系统的电子模型,观测并记录在不同阻尼比ζ时的跃响应曲线,并测定其超调量σ%及过渡过程时间ts。

三、实验原理1．一阶系统系统传递函数为：? s =r(s)=ts+1 模拟运算电路如图1- 1所示：c(s)k图 1- 1一阶系统模拟电路图由图 1-1得u0(s)r2/r1k==i2在实验当中始终取r2= r1，则k=1，t= r2c 取不同的时间常数t分别为： 0.25、 0.5、1.0。

记录不同时间常数下阶跃响应曲线，测量并记录其过渡过程时间ts，将参数及指标填在表1-1内。

2．二阶系统:系统传递函数为：? s =r(s)=s+2ζω1-2所示：c(s)ωn2ns+ωnωn=1弧度/秒，则系统结构如图图1-2二阶系统结构图根据结构图，建立的二阶系统模拟线路如图1-3所示：图1-3二阶系统模拟电路图取r2c1=1 ，r3c2 =1，则r4=r4c2=2ξ及ξ=2r3r114c2ζ取不同的值ζ=0.25 , ζ=0.5 , ζ=1.0，观察并记录阶跃响应曲线，测量超调量σ%，计算过渡过程时间ts，将参数及指标填入表1-2内。

四、实验设备1.hhmn-1电子模拟机一台。

2.pc机一台。

3.数字式万用表一块。

五、实验步骤1.熟悉hhmn-1型电子模拟机的使用方法，将各运算放大器接成比例器，通电调零。

2.断开电源，按照实验说明书上的条件和要求，计算电阻和电容的取值，按照模拟线路图搭接线路，不用的运算放大器接成比例器。