计算机仿真(龙格库塔方法)的软件VB设计与实现
matlab龙格库塔方法求解二元二阶常微分方程组
matlab龙格库塔方法求解二元二阶常微分方程组文章标题：深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用摘要：在科学与工程领域，常常需要求解复杂的微分方程组，而matlab作为一种强大的数学工具，提供了许多求解微分方程组的方法。
本文将深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用，以便读者全面理解该方法并能灵活应用于实际问题中。
正文：一、介绍龙格库塔方法龙格-库塔法（Runge-Kutta methods）是一种数值求解常微分方程的方法，通过将微分方程的解进行离散化，将微分方程转化为差分方程，从而进行数值求解。
龙格库塔方法通过迭代计算，能够得到微分方程的数值解，广泛应用于科学计算和工程技术领域。
二、matlab中的龙格库塔方法在matlab中，龙格库塔方法通过ode45函数实现，该函数能够对一阶或高阶常微分方程进行数值求解。
用户可以通过设定初始条件、微分方程表达式，以及积分区间等参数，快速得到微分方程的数值解。
ode45函数采用自适应步长的方式进行求解，能够有效解决微分方程解的数值稳定性和精确度问题。
三、龙格库塔方法在求解二元二阶常微分方程组中的应用考虑如下形式的二元二阶常微分方程组：$$\begin{cases}y_1' = f_1(t, y_1, y_2) \\y_2' = f_2(t, y_1, y_2)\end{cases}$$其中，$y_1(t)$和$y_2(t)$是未知函数，$f_1(t, y_1, y_2)$和$f_2(t,y_1, y_2)$分别表示其对应的函数表达式。
通过matlab中的ode45函数，可以将该二元二阶常微分方程组转化为一阶常微分方程组的形式，然后利用龙格库塔方法进行数值求解。
设定初始条件$y_1(0) = y1_0, y_2(0) = y2_0$，对应的一阶方程组为：$$\begin{cases}u_1' = u_3 \\u_2' = u_4 \\u_3' = f_1(t, u_1, u_2) \\u_4' = f_2(t, u_1, u_2)\end{cases}$$其中，$u_1(t) = y_1(t), u_2(t) = y_2(t), u_3(t) = y_1'(t), u_4(t) =y_2'(t)$，通过ode45函数求解该一阶常微分方程组即可得到原二元二阶常微分方程组的数值解。
数控计算机仿真课程设计最终版(VB)
广东工业大学课程设计任务书题目名称 数控系统的计算机仿真实现学生学院 机电工程学院 专业班级 机械设计制造及其自动化姓 名（学号）一、课程设计的内容对于给定的一段NC 代码，用VB 或其他高级语言编写程序解释、插补，在PC 机上仿真数控装置，进行图形描绘、坐标值显示、步进电机控制模拟显示及信号输出、冷却液和主轴开关量控制模拟显示及信号输出。
二、课程设计的要求与数据具体要求如下：(1) NC 代码中包含的代码类型有：G90 G54（G92） G00 G01 G02 G03 M03 M05 M08 M09 M30 例：下面给出一个具体的图形示意图，NC 代码及其加工轨迹图：% O0000 N106G0G90G54X10.Y20.M03M08N108Z50. N110Z10.N112G1Z-1. N114Y15.0 N118G2X15.Y13.09J7.5 N120X20.Y15.I-5.0J5.59 N122G1Y20.0 N126X10. N128G0Z50. N130M5M09 N136M30 %(2)、要求根据NC 代码屏幕模拟加工过程，图形显示位置，坐标值显示，辅助功能状态显示（冷却液和主轴开关量控制模拟显示）。
(3)、PC 机模拟加工过程中，要求有实时的驱动三轴步进电机的控制信号、控制冷却液和主轴转动的开关图1 工件平面图量输出控制信号。
假设信号从计算机并行打印口的数据信号线输出，端口地址为0x378。
并行口数据线分配如下（低电平有效）：表一并行口数据线信号定义数据线信号D0 D1 D2 D3 D4 D5 D6 D7定义PulseX DirX PulseY DirY PulseZ DirZ 主轴控制信号冷却液控制信号三、课程设计应完成的工作每个学生应在规定时间内，独立完成所选题目。
运用VB编程语言，编写计算机软件在WINDOWS实现数控装置的计算机仿真。
要求清楚地分析问题、提出算法、确定人机界面、列出流程图，最后用程序验证，完成软件测试，并且提交程序说明书。
计算机仿真习题及答案
计算机仿真试题1．编写一个函数，使其能够产生如下的分段函数：错误！未找到引用源。
并调用此函数，绘制x=[0,+2]范围内的f(x)*f(x+2) 。
（10分）function y=f(x)if x<=2y=0.5*x;else if x>6y=0.5;else y=1.5-0.25*x;endendx=0:0.05:2;y= f(x)’*f(x+2));plot(x,y)图 1-12．已知4阶龙格-库塔算法如下：试利用该算法求解以下微分方程：（15分）本题可以调用MATLAB函数中龙格-库塔算法函数ode45，首先编写m文件：function dy=func(x,y)dy=-y+1;end再在主窗口调用此文件：[x,y]=ode45('func',[0,5],0)%这里的[0,5]为任取区间，表示方程在此范围的解。
运行结果如下:x =0.00010.00010.00020.00020.00050.00070.00100.00120.00250.00370.00500.00620.01250.01880.02510.0313 0.06270.09410.12550.15690.28190.40690.53190.65690.78190.90691.03191.15691.28191.40691.53191.65691.78191.90692.03192.15692.28192.40692.53192.65692.78192.90693.03193.15693.28193.40693.53193.65693.78193.90694.03194.15694.28194.40694.53194.65694.74274.82854.91425.0000y =0.00010.00010.00020.00020.00050.00070.00100.00120.0025 0.0037 0.0050 0.0062 0.0124 0.0186 0.0248 0.0309 0.0608 0.0898 0.1180 0.1452 0.2457 0.33430.41250.48160.54250.59630.64370.68550.72250.75510.78390.80930.83170.85150.86890.88430.89790.90990.92050.92980.93810.94540.95180.95740.96240.96690.97080.97420.97720.97990.98230.98430.98620.98780.98920.99050.99130.99200.99270.9933为只管起见，我们使用函数命令画出x-y（plot(x,y)）的关系如下图：图1-23．用matlab语言求下列系统的状态方程、传递函数、零极点增益、和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：（15分）（1）G（s）=324327242410355024s s ss s s s+++++++（2）.X=2.25 -5 -1.25 -0.542.25 -4.25 -1.25 -0.2520.25 -0.5 -1.25 -121.25 -1.75 -0.25 -0.75 0X⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦uy= [0 2 0 2] X解：（1）a)求对应状态方程参数：num=[1 07 24 24]; den=[1 10 35 50 24]; [A,B,C,D]=tf2ss(num,den) 运行结果：A =-10 -35 -50 -241 0 0 00 1 0 00 0 1 0B =1C =1 7 24 24D =故，状态方程为：.X = x+ uY=[1 7 24 24]xb)求对应零极点增益模型参数：num=[1 07 24 24]; den=[1 10 35 50 24]; [Z,P,K]=tf2zp(num,den) 运行结果如下： Z =-2.7306 + 2.8531i -2.7306 - 2.8531i -1.5388P = -4.0000 -3.0000 -2.0000 -1.0000K = 1故变换后的零极点模型为： G(s)=c)求对应部分分式型：num=[1 07 24 24]; den=[1 10 35 50 24]; [R,P,H]=residue(num,den) 运行结果如下： R =4.0000 -6.0000 2.0000 1.0000P =-4.0000 -3.0000 -2.0000 -1.0000H = []故变换后的部分分式模型为：11223644)(+++++-+=s s s s s G（2）由题给条件，知：A=[2.25 -5 -1.25 -0.5; 2.25 -4.25 -1.25 -0.25;0.25 -0.5 -1.25 -1;1.25 -1.75-10 -35 -50 -24 1 0 0 0 0 1 0 0 0 0 1 010 0 0-0.25 -0.75] B=[4;2;2;0] C=[0 2 0 2],D=0 a)求传递函数矩阵： [num,den]=ss2tf(A,B,C,D) 运行结果为： num =0 4.0000 14.0000 22.0000 15.0000 den =1.0000 4.0000 6.2500 5.25002.2500 故，所对应传递函数模型为：25.225.525.641522144)(23423+++++++=s s s s s s s s Gb)求零极点模型：num=[0 4 14 22 15];en=[1 4 6.25 5.25 2.25]; [Z,P,K]=tf2zp(num,den) 运行结果为： Z =-1.0000 + 1.2247i -1.0000 - 1.2247i -1.5000 P =-1.5000 -1.5000 -0.5000 + 0.8660i -0.5000 - 0.8660iK =4.0000故，零极点模型为：)866.05.0()5.1()2247.11)(5.1(4)(2i s s i s s s G ±++±++=c)求对应部分分式模型： [R,P,H]=residue(num,den) 运行结果为： R =4.0000 -0.0000-0.0000 - 2.3094i -0.0000 + 2.3094iP =-1.5000 -1.5000 -0.5000 + 0.8660i -0.5000 - 0.8660iH = []故变换后的部分分式模型为：i s ii s i s s G 866.05.03094.2866.05.03094.25.14)(+++-+-++=4．已知一单位反馈系统开环传递函数为：，试绘制系统Nyquist图，判断闭环系统的稳定性，并求其单位阶跃响应。
龙格库塔方法及其matlab实现
龙格库塔方法及其matlab实现龙格-库塔方法及其matlab实现摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具达到求解目的。
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。
MatLab软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。
MatLab是英文MATrix LABoratory （矩阵实验室）的缩写。
在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。
关键词：龙格-库塔 matlab 微分方程1.前言1.1：知识背景龙格－库塔法（Runge-Kutta）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
通常所说的龙格库塔方法是相对四阶龙格库塔而言的，成为经典四阶龙格库塔法。
该方法具有精度高，收敛，稳定，计算过程中可以改变步长不需要计算高阶导数等优点，但是仍需计算在一些点上的值，比如四阶龙格-库塔法没计算一步需要计算四步，在实际运用中是有一定复杂性的。
Matlab是在20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB。
经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。
从这时起，MATLAB的内核采用C语言编写，而且除原有的数值计算能力外，还新增了数据图视功能。
MATLAB以商品形式出现后，仅短短几年，就以其良好的开放性和运行的可靠性，使原先控制领域里的封闭式软件包（如英国的UMIST，瑞典的LUND和SIMNON，德国的KEDDC）纷纷淘汰，而改以MATLAB为平台加以重建。
基于vb与matlab的数学与数字信号处理实验系统的设计与实现
基于VB与MATLAB的数学与数字信号处理实验系统的设计与实现一、引言随着科技的不断发展，数学与数字信号处理技术在各个领域的应用越来越广泛。
为了更好地满足教学和科研的需求，我们设计并实现了一个基于VB（Visual Basic）与MATLAB的数学与数字信号处理实验系统。
该系统能够提供丰富的实验环境，帮助学生和研究者进行数学与数字信号处理的学习、研究和实验。
二、系统设计1. 系统架构：本系统采用VB作为前端开发工具，实现用户界面和实验管理功能；后端使用MATLAB进行数学与数字信号处理算法的实现。
通过VB与MATLAB 的接口，实现前后端的交互。
2. 实验管理模块：该模块包括实验的创建、编辑、删除和运行等功能。
用户可以根据需要创建新的实验，对实验进行编辑和删除，以及运行实验并查看结果。
3. 算法实现模块：该模块利用MATLAB强大的数学与数字信号处理功能，实现各种算法。
包括但不限于：信号生成、滤波器设计、频谱分析、信号调制解调等。
4. 用户界面模块：该模块使用VB设计，提供直观、易用的用户界面。
用户可以通过界面进行实验管理，选择算法，输入参数，查看结果等操作。
三、系统实现1. VB与MATLAB接口：通过VB与MATLAB的接口，实现VB与MATLAB之间的数据交换。
VB将用户输入的参数传递给MATLAB，MATLAB进行算法计算后将结果返回给VB，并在界面上显示。
2. 算法实现：利用MATLAB的函数库，实现各种数学与数字信号处理算法。
对于每个算法，我们编写相应的MATLAB函数，并通过VB调用这些函数来执行实验。
3. 用户界面设计：使用VB设计友好的用户界面，包括实验管理界面、算法选择界面、参数设置界面和结果显示界面等。
用户可以通过界面方便地进行实验操作和管理。
四、结论本系统是一个功能强大、易于使用的数学与数字信号处理实验平台。
通过VB与MATLAB的结合，我们实现了实验管理的便捷性和算法实现的强大功能。
2024年度-《VB程序设计教程》PPT课件(全)
If...Then...Else语句
03
14
控制结构
Select Case语句
1
循环结构
2
For...Next循环
3
15
控制结构
Do...Loop循环 While...Wend循环
16Leabharlann 03VB界面设计17
窗体设计
窗体的类型与属性
介绍VB中不同类型的窗体，如标准窗体、MDI窗体等，以及窗 体的常用属性，如名称、标题、大小、位置等。
跨平台支持
支持跨平台开发，可运 行在Windows、Linux、Mac
等操作系统上。
5
VB的应用领域
桌面应用程序开发
利用VB可以快速开发Windows桌面应用程 序。
数据库应用开发
VB提供强大的数据库访问功能，可用于开发 数据库应用程序。
Web应用程序开发
通过可以开发 Web应用程 序。
运行时错误
根据错误提示信息定位问题所在，检查相关变量的值、函数调用等 ，修复逻辑错误或资源访问问题。
逻辑错误
仔细分析算法逻辑和程序流程，通过添加调试输出或日志记录来帮 助理解问题所在，并进行相应的逻辑调整。
35
程序优化策略探讨
代码优化
消除冗余代码、减少不必要的变量和对象创 建、优化循环结构等，提高代码执行效率。
多媒体应用开发
利用VB可以开发多媒体应用程序，如音频、 视频处理等。
6
02
VB编程基础
7
数据类型与变量
数值型
Integer、Long、Single、Double等
字符串型
String
8
数据类型与变量
布尔型
Boolean
四阶龙格-库塔(R-K)方法求常微分方程
中南大学MATLAB程序设计实践材料科学与工程学院2013年3月26日一、编程实现“四阶龙格－库塔（R-K）方法求常微分方程”，并举一例应用之。
【实例】采用龙格-库塔法求微分方程：,0)(1'00x x y y y 1、算法说明：在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。
其算法公式为：)22(63211k k k h y y nn其中：),()21,21()21,21(),(3423121hk y h x f k hk y h x f k hk y h x f k y x f k n nn n n n n n 2、流程图：2.1、四阶龙格－库塔（R-K ）方法流程图：2.2、实例求解流程图：输入待求微分方程、求解的自变量范围、初值以及求解范围内的取点数等。
确定求解范围内的步长k = 取点数？否求解：),()21,21()21,21(),(3423121hk y h x f k hk y h x f k hk y h x f k y x f k n nnn n n n n 求解并输出：)22(63211k k k h y y nn是结束算法开始输入求解的自变量范围求出待求简单微分方程的真值解用MA TLAB自带函数ode23求解待求微分方程用自编函数四阶龙格－库塔（R-K）方法求解待求微分方程结束3、源程序代码3.1、四阶龙格－库塔（R-K）方法源程序：function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)%Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x))%此程序可解高阶的微分方程。
只要将其形式写为上述微分方程的向量形式%函数 f(x,y): fun%自变量的初值和终值：x0, xt%y0表示函数在x0处的值，输入初值为列向量形式%自变量在[x0,xt]上取的点数：PointNum%varargin为可输入项，可传适当参数给函数f(x,y)%x：所取的点的x值%y：对应点上的函数值if nargin<4 | PointNum<=0PointNum=100;endif nargin<3y0=0;endy(1,:)=y0(:)'; %初值存为行向量形式h=(xt-x0)/(PointNum-1); %计算步长x=x0+[0:(PointNum-1)]'*h; %得x向量值for k=1:(PointNum) %迭代计算f1=h*feval(fun,x(k),y(k,:),varargin{:});f1=f1(:)'; %得公式k1 f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:});f2=f2(:)'; %得公式k2 f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:});f3=f3(:)'; %得公式k3 f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:});f4=f4(:)'; %得公式k4 y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得y(n+1) end3.2、实例求解源程序：%运行四阶R-K法clear, clc %清除内存中的变量x0=0;xt=2;Num=100;h=(xt-x0)/(Num-1);x=x0+[0:Num]*h;a=1;yt=1-exp(-a*x); %真值解fun=inline('-y+1','x','y'); %用inline构造函数f(x,y)y0=0; %设定函数初值PointNum=5; %设定取点数[x1,y1]=ode23(fun,[0,2],0);[xr,yr]=MyRunge_Kutta(fun,x0,xt,y0,PointNum);MyRunge_Kutta_x=xr'MyRunge_Kutta_y=yr'plot(x,yt,'k',x1,y1,'b--',xr,yr,'r-')legend('真值','ode23','Rung-Kutta法解',2)hold onplot(x1,y1,'bo',xr,yr,'r*')4、程序运行结果：MyRunge_Kutta_x =0 0.5000 1.0000 1.5000 2.0000MyRunge_Kutta_y =0 0.3932 0.6318 0.7766 0.8645二、变成解决以下科学计算问题：（一）[例7-2-4] 材料力学复杂应力状态的分析——Moore 圆。
Matlab中龙格-库塔(Runge-Kutta)方法原理及实现
函数功能编辑本段回目录ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。
ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)³。
解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解.使用方法编辑本段回目录[T,Y] = ode45(odefun,tspan,y0)odefun 是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名tspan 是区间[t0 tf] 或者一系列散点[t0,t1,...,tf]y0 是初始值向量T 返回列向量的时间点Y 返回对应T的求解列向量[T,Y] = ode45(odefun,tspan,y0,options)options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等[T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options)在设置了事件参数后的对应输出TE 事件发生时间YE 事件解决时间IE 事件消失时间sol =ode45(odefun,[t0 tf],y0...)sol 结构体输出结果应用举例编辑本段回目录1 求解一阶常微分方程程序:一阶常微分方程odefun=@(t,y) (y+3*t)/t^2; %定义函数tspan=[1 4]; %求解区间y0=-2; %初值[t,y]=ode45(odefun,tspan,y0);plot(t,y) %作图title('t^2y''=y+3t,y(1)=-2,1<t<4')legend('t^2y''=y+3t')xlabel('t')ylabel('y')% 精确解% dsolve('t^2*Dy=y+3*t','y(1)=-2')% ans =一阶求解结果图% (3*Ei(1) - 2*exp(1))/exp(1/t) - (3*Ei(1/t))/exp(1/t)2 求解高阶常微分方程关键是将高阶转为一阶,odefun的书写.F(y,y',y''...y(n-1),t)=0用变量替换,y1=y,y2=y'...注意odefun方程定义为列向量dxdy=[y(1),y(2)....]程序:function Testode45tspan=[3.9 4.0]; %求解区间y0=[2 8]; %初值[t,x]=ode45(@odefun,tspan,y0);plot(t,x(:,1),'-o',t,x(:,2),'-*')legend('y1','y2')title('y'' ''=-t*y + e^t*y'' +3sin2t')xlabel('t')ylabel('y')function y=odefun(t,x)y=zeros(2,1); % 列向量y(1)=x(2);y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);endend高阶求解结果图相关函数编辑本段回目录ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tbMatlab中龙格-库塔(Runge-Kutta)方法原理及实现(自己写的，非直接调用)龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。
基于VB的常用水力学公式的程序设计34
基于VB的常用水力学公式的程序设计摘要：在水力学中，有很多公式的求解涉及到非线性方程的求解问题。
采用人工手算不仅工作量大，且精度不高。
本文运用Visual Basic编程语言，采用试算迭代法、四阶龙格库塔法完成了对2个常用水力学公式，即：正常水深计算公式及水平不透水层上的均质土坝单宽渗流量和逸出点高度计算公式的程序设计，并且编写了相应的应用程序界面，方便运用计算机计算,提高了这2个常用水力学公式计算的速度和精度。
关键词：水力学；公式；迭代法；四阶龙格库塔法；VB；程序设计1 Visual Basic简介Visual Basic(简称“VB”)是microsoft公司1991年推出功能很强面向对象的语言，具有友好图形用户界面的可视化高级计算机程序开发语言。
vb是准编译语言，它产生的最终代码不是独立可执行，而是需要一个动态链接库区解释才能执行。
vb进行程序开发步骤是：首先设计界面，然后编写代码。
设计界面时，可以将需要的控件按照一定的位置放进窗口中，而编写代码是让计算机完成给定的任务。
2 利用vb编程计算水力学公式明渠中由于渠道横断面的几何形状或尺寸，粗糙度或底坡沿程改变，或在明渠中修建人工筑物（闸、桥梁、涵洞）等都会改变水流的均匀状态，造成水深和流速等水力要素沿程改变，从而产生非均匀流动。
人工渠道或天然河道中的水流大多数都属于非均匀流，因此研究明渠恒定非均匀流对解决生产实际问题有重要意义。
2.1正常水深计算由均匀流关系，推导得（2.1.1）对梯形断面，过水断面面积与湿周为；（2.1.2）当取边坡系数m为零时，为矩形断面相应的水力要素。
将梯形渠道面积与湿周公式代入（2.1.1）式可导出已知梯形渠道底宽b求正常水深h的m=0迭代公式如下：（2.1.3）迭代计算水深时可取初值为将式（2.1.3）中边坡系数以m=0替换时，可得出矩形渠道正常水深h迭代公式为：以上各式中j表示迭代次数。
具体的迭代过程为，先计算出式中的常数项，然后任设一初值h(0)，代入迭代式右边的3结语本文运用Visual Basic程序设计语言，以迭代试算法和四阶龙格库塔法为基本计算方法，完成了2个常用水力学公式，即：正常水深计算公式及水平不透水层上的均质土坝单宽渗流量q和逸出点高度a0的公式的程序设计。
哈工大 计算机仿真技术实验报告 实验3 利用数值积分算法的仿真实验
模型的稳定性：当步距 h=5.0e-5 时，前向欧拉法和后向欧拉法明显失真， 随着步距的减小， 二阶显式 Adams 法， 梯形法和显式四阶 Runge-Kutta 法的波形 变化不大，而前向欧拉法和后向欧拉法的波形得到明显改善。所以显式四阶 Runge-Kutta 法，二阶显式 Adams 法和梯形法的稳定性较好，前向欧拉法和后向 欧拉法的稳定性较差。 模型的精度和离散时间间隔：步距为 h=5.0e-6 时，显式四阶 Runge-Kutta 法 精度最高，其次是二阶显式 Adams 法和梯形法。步距为 h=5.0e-7 时，前向欧拉 法和后向欧拉法仿真精度才达到要求。所以，显式四阶 Runge-Kutta 法，二阶显 式 Adams 法和梯形法模型的精度较高，离散时间间隔要求低，其中，显式四阶 Runge-Kutta 法模型的精度最高，其次是二阶显式 Adams 法，由于是二次函数较 复杂，函数曲线与真实曲线较为接近；再次精确的是梯形法，取梯形面积，误差 也较小；前向欧拉法和后向欧拉法模型的精度较低，由于取的是矩形面积，离散 时间间隔要求高。
实验 3 利用数值积分算法的仿真实验
(
一、 实验目的
1) 熟悉 MATLAB 的工作环境；
2) 掌握 MATLAB 的 .M 文件编写规则，并在命令窗口调试和运行程序； 3) 掌握利用欧拉法、梯形法、二阶显式 Adams 法及四阶龙格库塔法构建系 统仿真模型的方法,并对仿真结果进行分析。
二、实验内容
上对应的标题。
四、实验原理
在连续系统的数字仿真算法中,较常用的有欧拉法、 梯形法、 二阶显式 Adams 法及显式四阶 Runge-Kutta 法等。欧拉法、梯形法和二阶显式 Adams 法是利用离 散相似原理构造的仿真算法，而显式四阶 Runge-Kutta 法是利用 Taylor 级数匹配 原理构造的仿真算法。 对于线性系统，其状态方程表达式为：
二维动力学微分方程迭代解法matlab龙格库塔
1. 二维动力学微分方程的重要性二维动力学微分方程是描述动力学系统中两个变量之间相互作用的数学模型。
这些方程在物理学、工程学、生物学等领域有着广泛的应用。
解决这些方程对于深入理解动力学系统的行为以及预测其未来的状态具有至关重要的意义。
2. 迭代解法的基本原理迭代解法是一种通过不断逼近真实解的方法，其基本原理是从一个初始猜测值开始，通过不断的重复计算来逼近真实解。
这种方法在解决复杂的微分方程时具有很高的实用性和灵活性。
3. MATLAB在迭代解法中的应用MATLAB是一种强大的数学计算软件，它提供了丰富的工具和函数来解决各种数学问题。
在迭代解法中，MATLAB提供了龙格-库塔（Runge-Kutta）方法来解决微分方程，可以高效地求解二维动力学微分方程。
4. 龙格-库塔方法的特点龙格-库塔方法是一种经典的迭代解法算法，它通过多步迭代来逼近微分方程的真实解。
这种方法在数值计算领域被广泛应用，其精度和稳定性都很高，并且适用于各种复杂的微分方程。
5. 在MATLAB中使用龙格-库塔方法解决二维动力学微分方程的步骤a. 定义微分方程的函数表达式b. 设置初始条件c. 调用MATLAB的龙格-库塔函数进行迭代计算d. 获取最终的解并进行结果分析6. 实例分析为了更好地理解在MATLAB中使用龙格-库塔方法解决二维动力学微分方程的过程，我们以一个具体的实例来进行分析。
7. 结论通过使用MATLAB中的龙格-库塔方法，我们可以高效地解决各种复杂的二维动力学微分方程，得到准确的数值解，从而更好地理解动力学系统的行为。
这为动力学系统的建模和分析提供了重要的工具和方法。
8. 实例分析假设我们有一个二维动力学系统，其微分方程表示为：\[ \begin{cases} \frac{dx}{dt} = y \\ \frac{dy}{dt} = -x - 0.5y\end{cases} \]我们希望使用MATLAB的龙格-库塔方法来求解该微分方程，并分析系统的行为。
计算机仿真（龙格库塔方法）的软件VB设计与实现
计算机仿真（龙格库塔⽅法）的软件VB设计与实现计算机仿真（龙格库塔⽅法）的软件VB设计与实现Dim a(0 To 10) As Single, b(0 To 10) As Single, c(0 To 10) As Single, d(0 To 10) As SingleDim e(0 To 10) As Single, h(0 To 10) As Single, p(0 To 10) As Single, q(0 To 10) As SingleDim n1 As Byte, n2 As Byte 'n1表⽰的是y 的阶数,n2表⽰的是输⼊函数的阶数Dim i As IntegerDim f(0 To 2000) As Single 'f(m)=y,也就是各个时刻的y值Dim X0(0 To 12) As Single, X1(0 To 12) As SingleDim dt As Single, u As Single 'dt为采样周期,u为输⼊Dim q1 As Single, p1 As Single, h1 As Single, e1 As Single Dim b0(0 To 10) As Single, bn(-4 To 10) As SingleDim m As IntegerDim max As SinglePrivate Sub Combo1_Click()n1 = Combo1.ListIndex + 1 'n1表⽰的是y的阶数For i = n1 + 1 To 8Text1(i).Visible = FalseLabel1(i).Visible = FalseLabel9(i).Visible = FalseNextFor i = 0 To n1Text1(i).Visible = TrueLabel1(i).Visible = TrueLabel9(i).Visible = TrueNextCombo2.ClearCombo2.Text = "请选择输出u的最⼤阶数"For i = 0 To n1Combo2.AddItem ((i) & "阶")NextEnd SubPrivate Sub Combo2_Click()n2 = Combo2.ListIndex 'n2表⽰的是输⼊函数的阶数For i = n2 + 1 To 8Text2(i).Visible = FalseLabel2(i).Visible = FalseLabel10(i).Visible = FalseNextFor i = 0 To n2Text2(i).Visible = TrueLabel2(i).Visible = TrueLabel10(i).Visible = TrueNextEnd SubPrivate Sub Calculate_Click()dt = Val(Text3.Text) 'dt是采样周期u = Val(Text4.Text)n1 = Combo1.ListIndex + 1If n1 = 0 Then '对付忘记选择最⾼阶数时的情况Call MsgBox("请选择阶数！", 48, "未选择阶数") Exit SubEnd IfFor i = 0 To 12 '设X0的初值都是0，且i⼤于8，这⼀点很重要！X0(i) = 0NextFor i = 0 To n1 - 1a(i) = -Val(Text1(i).Text) / Val(Text1(n1).Text) 'a(i)b(i)c(i)d(i)等等都是计算的中间变量NextFor i = 0 To n2b0(i) = Val(Text2(i).Text) / Val(Text1(n1).Text)NextIf Combo3.ListIndex = 1 Thenn1 = n1 + 1For i = n1 - 1 To 1 Step -1a(i) = a(i - 1)Nexta(0) = 0End IfFor i = -4 To -1bn(i) = 0Nextb(0) = a(0) * a(n1 - 1)For i = 1 To n1 - 1b(i) = a(i - 1) + a(i) * a(n1 - 1)Nextc(0) = a(0) * b(n1 - 1)For i = 1 To n1 - 1c(i) = b(i - 1) + a(i) * b(n1 - 1)Nextd(0) = a(0) * c(n1 - 1)For i = 1 To n1 - 1d(i) = c(i - 1) + a(i) * c(n1 - 1)NextFor i = 0 To n1 - 1e(i) = dt * a(i) + dt * dt / 2 * b(i) + dt * dt * dt / 6 * c(i) + dt * dt * dt * dt / 24 * d(i)h(i) = dt * dt / 2 * a(i) + dt * dt * dt / 6 * b(i) + dt * dt * dt * dt / 24 * c(i)p(i) = dt * dt * dt / 6 * a(i) + dt * dt * dt * dt / 24 * b(i)q(i) = dt * dt * dt * dt / 24 * a(i)NextIf n2 = 0 Then 'n2为0和不为0的公式是不⼀样的，表⽰输⼊函数u为n2阶'b0(0) = Val(Text2(0).Text) / Val(Text1(n1).Text)m = 1Doq1 = 0: p1 = 0: h1 = 0: e1 = 0For i = 0 To n1 - 1 '这些个也是中间变量q1 = q1 + q(i) * X0(i + 1)p1 = p1 + p(i) * X0(i + 1)h1 = h1 + h(i) * X0(i + 1)e1 = e1 + e(i) * X0(i + 1)NextFor i = 1 To n1X1(i) = X0(i) + dt * X0(i + 1) + dt * dt / 2 * X0(i + 2) + dt * dt * dt / 6 * X0(i + 3) + dt * dt * dt * dt / 24 * X0(i + 4)NextX1(n1) = X1(n1) + e1 + (dt + dt * dt / 2 * a(n1 - 1) + dt * dt * dt / 6 * b(n1 - 1) + dt * dt * dt * dt / 24 * c(n1 - 1)) * b0(0) * u If n1 >= 4 ThenX1(n1 - 3) = X1(n1 - 3) + q1 + dt * dt * dt * dt / 24 * b0(0) * uX1(n1 - 2) = X1(n1 - 2) + p1 + (dt * dt * dt / 6 + dt * dt * dt * dt / 24 * a(n1 - 1)) * b0(0) * uX1(n1 - 1) = X1(n1 - 1) + h1 + (dt * dt / 2 + dt * dt * dt / 6 * a(n1 - 1) + dt * dt * dt * dt / 24 * b(n1 - 1)) * b0(0) * uElseIf n1 >= 3 ThenX1(n1 - 2) = X1(n1 - 2) + p1 + (dt * dt * dt / 6 + dt * dt * dt * dt / 24 * a(n1 - 1)) * b0(0) * uX1(n1 - 1) = X1(n1 - 1) + h1 + (dt * dt / 2 + dt * dt * dt / 6 * a(n1 - 1) + dt * dt * dt * dt / 24 * b(n1 - 1)) * b0(0) * uElseIf n1 >= 2 ThenX1(n1 - 1) = X1(n1 - 1) + h1 + (dt * dt / 2 + dt * dt * dt / 6 * a(n1 - 1) + dt * dt * dt * dt / 24 * b(n1 - 1)) * b0(0) * uEnd Iff(m) = X1(1)For i = 1 To n1X0(i) = X1(i)Nextm = m + 1Loop Until m = 2000 '循环计算2000次ElseFor i = (n2 + 1) To 10b0(i) = 0NextFor i = 0 To n1 - 1bn(i) = b0(i) + b0(n1) * a(i)Nextm = 1DoFor i = 1 To n1X1(i) = X0(i) + dt * X0(i + 1) + dt * dt / 2 * X0(i + 2) + dt * dt * dt / 6 * X0(i + 3) + dt * dt * dt * dt / 24 * X0(i + 4) + (h(n1 - i) * bn(n1 -1) + p(n1 - i) * bn(n1 - 2) + q(n1 - i) * bn(n1- 3)) * uNextFor i = 1 To n1X1(i) = X1(i) + e(n1 - i) * X0(1) + h(n1 - i) * X0(2) + p(n1 - i) * X0(3) + q(n1 - i) * X0(4) + (dt * bn(n1 - i) + dt * dt / 2 * bn(n1 - 1 - i) + dt * dt * dt / 6 * bn(n1 - 2 - i) + dt * dt * dt * dt / 24 * bn(n1 - 3 - i)) * uNextf(m) = X1(1)For i = 1 To n1X0(i) = X1(i)Nextm = m + 1Loop Until m = 2000End Ifmax = f(1)For m = 2 To 1999If Abs(f(m)) > Abs(max) Then '开始绘图max = f(m)End IfNext'If max > 0 ThenPicture1.ClsPicture1.Scale (-20 * dt, 1.2 * max)-(2020 * dt, -0.05 * max) Picture1.Line (0, 0)-(2002 * dt, 0)Picture1.Line (0, 0)-(0, 1.18 * max)Picture1.Line (2002 * dt, 0)-(1960 * dt, -0.025 * max) Picture1.Line (2002 * dt, 0)-(1960 * dt, 0.025 * max)Picture1.Line (0, 1.18 * max)-(10 * dt, 1.1 * max)Picture1.Line (0, 1.18 * max)-(-10 * dt, 1.1 * max)'Else' Picture1.Cls' Picture1.Scale (-5 * dt, 0.2 * Abs(max))-(2002 * dt, 1.2 * max) ' Picture1.Line (0, 0)-(2002 * dt, 0)' Picture1.Line (0, 0)-(0, 1.2 * max)' Picture1.Line (2002 * dt, 0)-(1900 * dt, 0.1 * max)' Picture1.Line (0, 1.2 * max)-(50 * dt, 1.1 * max)' End IfFor m = 2 To 1999Picture1.Line ((m - 1) * dt, f(m - 1))-(m * dt, f(m)), vbRed Next 'If (f(1999) - f(1800)) < 0.00001 Then'Picture1.Line (0, f(1999))-(2000 * dt, f(1999))' Picture1.CurrentX = 0'Picture1.CurrentY = f(1999) + 0.1 * max'Picture1.Print "y ∞=" & f(1999)'End IfIf Combo3.ListIndex = 0 ThenPicture1.Line (0, f(1999))-(2000 * dt, f(1999))Picture1.CurrentX = 0Picture1.CurrentY = f(1999) + 0.1 * maxPicture1.Print "y ∞=" & f(1999)ElsePicture1.Line (0, 0)-(2000 * dt, u * 2000 * dt)End IfEnd SubPrivate Sub Combo3_Click()Text4.Visible = TrueLabel12.Visible = TrueIf Combo3.ListIndex = 0 ThenLabel5.Caption = "请输⼊阶跃信号u的系数："ElseLabel5.Caption = "请输⼊斜坡信号u的系数：" End If End SubPrivate Sub dtxnzb_Click()Dim y00 As SingleDim thigema As SingleDim ts As SingleDim ess As SingleIf Combo3.ListIndex = 0 Theny00 = f(1999)thigema = max - y00m = 2000Dom = m - 1Loop Until Abs(f(m) - y00) >= (0.05 * y00)ts = m * dtess = u - y00Frame1.Visible = TrueLabel6.Caption = "超调量σ:" & thigemaLabel7.Caption = "调节时间ts:" & tsLabel8.Caption = "稳态误差ess:" & essElseFrame1.Visible = TrueLabel6.Caption = "稳态误差ess:" & (u * 1999 * dt - f(1999)) / (u * 1999 * dt) & "%" End If End Sub。
四阶龙格—库塔法在电机数字仿真中的适用范围与最佳步长的确定
四阶龙格—库塔法在电机数字仿真中的适用范围与最佳步长
的确定
邬学义
【期刊名称】《计算机仿真》
【年(卷),期】1992(000)003
【摘要】本文就四阶龙格——库塔法在电机数字仿真中的有关问题,包括电机方程病态性的判断与仿真时最佳步长的确定等,进行了探讨。
实践证明,应用本文提出的判据和方法,可使电机的数字仿真得到快速,准确的效果。
【总页数】4页(P26-29)
【作者】邬学义
【作者单位】无
【正文语种】中文
【中图分类】TM3
【相关文献】
1.四阶龙格-库塔法在捷联惯导系统姿态解算中的应用 [J], 张春慧;吴简彤;何昆鹏;周雪梅
2.四阶龙格—库塔法的原理及其应用 [J], 冯建强;孙诗一;
3.四阶龙格—库塔法在电机数字仿真中的应用问题 [J], 孟传富;张明玉
4.四阶龙格_库塔法在火控解算中的应用 [J], 李丹
5.变步长四阶龙格库塔方法在水文水力计算中的应用 [J], 梁建;李宝春;祁涛;卢俊
因版权原因，仅展示原文概要，查看原文内容请购买。
控制系统计算机仿真及辅助设计实验报告
阶跃
num=[0.8,0,-20];
den=[1,0,-40,0];
sys=tf(num,den);
t=0:0.01:1;
step(sys,t)
实验图形
室温控制系统校正装置设计
已知某室温控制系统为单位负反馈，某开环传递函数为： ，试用Bode图设计法对系统进行滞后串联校正设计，使系统满足;
系统在斜坡信号作用下，系统的速度误差系数 ≥30
(2)比较这几种方法：
对于四阶龙格-库塔方法
真值
1
0.9048
0.8187
0.7408
0.6703
0.6065
0.5488
0.4966
0.4493
0.4066
0.3679
龙库
1
0.9048
0.8187
0.7408
0.6703
0.6065
0.5488
0.4966
0.4493
0.4066
0.3679
误差
step(sys,t)
单位脉冲响应图像
单位阶跃响应图像
实验二
2-2.用MATLAB语言求下列系统的状态方程、传递函数、零极点增益、和部分分式形式的模型参数，并分别写出其相应的数学模型表达式：
1.G（s）=
2. =
Y=[0 2 0 2] X
1.解：（1）状态方程模型参数:
编写MATLAB程序如下
>> num=[1 7 24 24];
（1）m文件程序为h=0.1;
disp('函数的数值解为'); %显示‘’中间的文字%
disp('y=');%同上%
y=1;
龙格库塔法-原理及程序实现
龙格库塔法一、基本原理：“龙格-库塔(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。
计算机仿真技术实验报告
《计算机仿真技术》实验指导书实验一 状态空间模型的仿真一、实验目的通过实验，学习4阶龙格-库塔法的基本思路和计算公式，加深理解4阶龙格-库塔法的原理和稳定域。
加深理解仿真的稳定性，仿真步长对仿真精度的影响。
二、实验内容1、线性定常系统[]1112223332010002001010060000600x x x x x u y x x x x -⎡⎤⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=-+=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦&&&；)(1000)0()0()0(321t u x x x =⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡2、非线性系统()()()()()()()()xt rx t ax t y t yt sx t bx t y t =-⎧⎨=-+⎩&& 其中：r=0.001, a=2⨯10-6, s=0.01, b=1⨯10-6, x(0)=12000, y(0)=600。
三、实验原理运用SIMULINK 仿真工具进行实验。
四、实验设备和仪器微型计算机、MATLAB 软件。
Sources(信号源)，Sink(显示输出),Continuous(线性连续系统)，Discrete（线性离散系统)，Function & Table （函数与表格)，Math(数学运算)， Discontinuities (非线性)，Demo （演示）五、实验方法运行MA TLAB ，在MA TLAB 窗口中按SimuLink 按钮，启动SimuLink 库浏览器，在浏览器窗口上选create a new modem 命令，得到一个空模型，从Library: SimuLink 窗口中找到需要的模块，将这些模块拖到空模型窗口中。
将空模型窗口中的排好，并按要求连接。
在保存好的模型窗口中，选Simulation\Paramters 命令设置各模块的参数和仿真参数。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算机仿真（龙格库塔方法）的软件VB设计与实现Dim a(0 To 10) As Single, b(0 To 10) As Single, c(0 To 10) As Single, d(0 To 10) As SingleDim e(0 To 10) As Single, h(0 To 10) As Single, p(0 To 10) As Single, q(0 To 10) As SingleDim n1 As Byte, n2 As Byte 'n1表示的是y 的阶数,n2表示的是输入函数的阶数Dim i As IntegerDim f(0 To 2000) As Single 'f(m)=y,也就是各个时刻的y值Dim X0(0 To 12) As Single, X1(0 To 12) As SingleDim dt As Single, u As Single 'dt为采样周期,u为输入Dim q1 As Single, p1 As Single, h1 As Single, e1 As SingleDim b0(0 To 10) As Single, bn(-4 To 10) As SingleDim m As IntegerDim max As SinglePrivate Sub Combo1_Click()n1 = Combo1.ListIndex + 1 'n1表示的是y的阶数For i = n1 + 1 To 8Text1(i).Visible = FalseLabel1(i).Visible = FalseLabel9(i).Visible = FalseNextFor i = 0 To n1Text1(i).Visible = TrueLabel1(i).Visible = TrueLabel9(i).Visible = TrueNextCombo2.ClearCombo2.Text = "请选择输出u的最大阶数"For i = 0 To n1Combo2.AddItem ((i) & "阶")NextEnd SubPrivate Sub Combo2_Click()n2 = Combo2.ListIndex 'n2表示的是输入函数的阶数For i = n2 + 1 To 8Text2(i).Visible = FalseLabel2(i).Visible = FalseLabel10(i).Visible = FalseNextFor i = 0 To n2Text2(i).Visible = TrueLabel2(i).Visible = TrueLabel10(i).Visible = TrueNextEnd SubPrivate Sub Calculate_Click()dt = Val(Text3.Text) 'dt是采样周期u = Val(Text4.Text)n1 = Combo1.ListIndex + 1If n1 = 0 Then '对付忘记选择最高阶数时的情况Call MsgBox("请选择阶数！", 48, "未选择阶数")Exit SubEnd IfFor i = 0 To 12 '设X0的初值都是0，且i大于8，这一点很重要！X0(i) = 0NextFor i = 0 To n1 - 1a(i) = -Val(Text1(i).Text) / Val(Text1(n1).Text) 'a(i)b(i)c(i)d(i)等等都是计算的中间变量NextFor i = 0 To n2b0(i) = Val(Text2(i).Text) / Val(Text1(n1).Text)NextIf Combo3.ListIndex = 1 Thenn1 = n1 + 1For i = n1 - 1 To 1 Step -1a(i) = a(i - 1)Nexta(0) = 0End IfFor i = -4 To -1bn(i) = 0Nextb(0) = a(0) * a(n1 - 1)For i = 1 To n1 - 1b(i) = a(i - 1) + a(i) * a(n1 - 1)Nextc(0) = a(0) * b(n1 - 1)For i = 1 To n1 - 1c(i) = b(i - 1) + a(i) * b(n1 - 1)Nextd(0) = a(0) * c(n1 - 1)For i = 1 To n1 - 1d(i) = c(i - 1) + a(i) * c(n1 - 1)NextFor i = 0 To n1 - 1e(i) = dt * a(i) + dt * dt / 2 * b(i) + dt * dt * dt / 6 * c(i) + dt * dt * dt * dt / 24 * d(i)h(i) = dt * dt / 2 * a(i) + dt * dt * dt / 6 * b(i) + dt * dt * dt * dt / 24 * c(i)p(i) = dt * dt * dt / 6 * a(i) + dt * dt * dt * dt / 24 * b(i)q(i) = dt * dt * dt * dt / 24 * a(i)NextIf n2 = 0 Then 'n2为0和不为0的公式是不一样的，表示输入函数u为n2阶'b0(0) = Val(Text2(0).Text) / Val(Text1(n1).Text)m = 1Doq1 = 0: p1 = 0: h1 = 0: e1 = 0For i = 0 To n1 - 1 '这些个也是中间变量q1 = q1 + q(i) * X0(i + 1)p1 = p1 + p(i) * X0(i + 1)h1 = h1 + h(i) * X0(i + 1)e1 = e1 + e(i) * X0(i + 1)NextFor i = 1 To n1X1(i) = X0(i) + dt * X0(i + 1) + dt * dt / 2 * X0(i + 2) + dt * dt * dt / 6 * X0(i + 3) + dt * dt * dt * dt / 24 * X0(i + 4)NextX1(n1) = X1(n1) + e1 + (dt + dt * dt / 2 * a(n1 - 1) + dt * dt * dt / 6 * b(n1 - 1) + dt * dt * dt * dt / 24 * c(n1 - 1)) * b0(0) * uIf n1 >= 4 ThenX1(n1 - 3) = X1(n1 - 3) + q1 + dt * dt * dt * dt / 24 * b0(0) * uX1(n1 - 2) = X1(n1 - 2) + p1 + (dt * dt * dt / 6 + dt * dt * dt * dt / 24 * a(n1 - 1)) * b0(0) * uX1(n1 - 1) = X1(n1 - 1) + h1 + (dt * dt / 2 + dt * dt * dt / 6 * a(n1 - 1) + dt * dt * dt * dt / 24 * b(n1 - 1)) * b0(0) * uElseIf n1 >= 3 ThenX1(n1 - 2) = X1(n1 - 2) + p1 + (dt * dt * dt / 6 + dt * dt * dt * dt / 24 * a(n1 - 1)) * b0(0) * uX1(n1 - 1) = X1(n1 - 1) + h1 + (dt * dt / 2 + dt * dt * dt / 6 * a(n1 - 1) + dt * dt * dt * dt / 24 * b(n1 - 1)) * b0(0) * uElseIf n1 >= 2 ThenX1(n1 - 1) = X1(n1 - 1) + h1 + (dt * dt / 2 + dt * dt * dt / 6 * a(n1 - 1) + dt * dt * dt * dt / 24 * b(n1 - 1)) * b0(0) * uEnd Iff(m) = X1(1)For i = 1 To n1X0(i) = X1(i)Nextm = m + 1Loop Until m = 2000 '循环计算2000次ElseFor i = (n2 + 1) To 10b0(i) = 0NextFor i = 0 To n1 - 1bn(i) = b0(i) + b0(n1) * a(i)Nextm = 1DoFor i = 1 To n1X1(i) = X0(i) + dt * X0(i + 1) + dt * dt / 2 * X0(i + 2) + dt * dt * dt / 6 * X0(i + 3) + dt * dt * dt * dt / 24 * X0(i + 4) + (h(n1 - i) * bn(n1 - 1) + p(n1 - i) * bn(n1 - 2) + q(n1 - i) * bn(n1- 3)) * uNextFor i = 1 To n1X1(i) = X1(i) + e(n1 - i) * X0(1) + h(n1 - i) * X0(2) + p(n1 - i) * X0(3) + q(n1 - i) * X0(4) + (dt * bn(n1 - i) + dt * dt / 2 * bn(n1 - 1 - i) + dt * dt * dt / 6 * bn(n1 - 2 - i) + dt * dt * dt * dt / 24 * bn(n1 - 3 - i)) * uNextf(m) = X1(1)For i = 1 To n1X0(i) = X1(i)Nextm = m + 1Loop Until m = 2000End Ifmax = f(1)For m = 2 To 1999If Abs(f(m)) > Abs(max) Then '开始绘图max = f(m)End IfNext'If max > 0 ThenPicture1.ClsPicture1.Scale (-20 * dt, 1.2 * max)-(2020 * dt, -0.05 * max)Picture1.Line (0, 0)-(2002 * dt, 0)Picture1.Line (0, 0)-(0, 1.18 * max)Picture1.Line (2002 * dt, 0)-(1960 * dt, -0.025 * max)Picture1.Line (2002 * dt, 0)-(1960 * dt, 0.025 * max)Picture1.Line (0, 1.18 * max)-(10 * dt, 1.1 * max)Picture1.Line (0, 1.18 * max)-(-10 * dt, 1.1 * max)'Else' Picture1.Cls' Picture1.Scale (-5 * dt, 0.2 * Abs(max))-(2002 * dt, 1.2 * max)' Picture1.Line (0, 0)-(2002 * dt, 0)' Picture1.Line (0, 0)-(0, 1.2 * max)' Picture1.Line (2002 * dt, 0)-(1900 * dt, 0.1 * max)' Picture1.Line (0, 1.2 * max)-(50 * dt, 1.1 * max)' End IfFor m = 2 To 1999Picture1.Line ((m - 1) * dt, f(m - 1))-(m * dt, f(m)), vbRed Next'If (f(1999) - f(1800)) < 0.00001 Then'Picture1.Line (0, f(1999))-(2000 * dt, f(1999))' Picture1.CurrentX = 0'Picture1.CurrentY = f(1999) + 0.1 * max'Picture1.Print "y ∞=" & f(1999)'End IfIf Combo3.ListIndex = 0 ThenPicture1.Line (0, f(1999))-(2000 * dt, f(1999))Picture1.CurrentX = 0Picture1.CurrentY = f(1999) + 0.1 * maxPicture1.Print "y ∞=" & f(1999)ElsePicture1.Line (0, 0)-(2000 * dt, u * 2000 * dt)End IfEnd SubPrivate Sub Combo3_Click()Text4.Visible = TrueLabel12.Visible = TrueIf Combo3.ListIndex = 0 ThenLabel5.Caption = "请输入阶跃信号u的系数："ElseLabel5.Caption = "请输入斜坡信号u的系数：" End IfEnd SubPrivate Sub dtxnzb_Click()Dim y00 As SingleDim thigema As SingleDim ts As SingleDim ess As SingleIf Combo3.ListIndex = 0 Theny00 = f(1999)thigema = max - y00m = 2000Dom = m - 1Loop Until Abs(f(m) - y00) >= (0.05 * y00)ts = m * dtess = u - y00Frame1.Visible = TrueLabel6.Caption = "超调量σ:" & thigemaLabel7.Caption = "调节时间ts:" & tsLabel8.Caption = "稳态误差ess:" & essElseFrame1.Visible = TrueLabel6.Caption = "稳态误差ess:" & (u * 1999 * dt - f(1999)) / (u * 1999 * dt) & "%" End IfEnd Sub。