计算机仿真技术实验报告-实验三
计算机仿真实验报告-实验三
一、实验内容：实验三 利用欧拉法、梯形数法和二阶显式Adams 法对RLC 串联电路的仿真1前向欧拉法状态方程：Du CX y Bu AX X m +=+=+•1 然后根据前向欧拉法（其中h 为步长）•++=m m m hX X X 1即可得到系统的差分方程2后向欧拉法根据前向欧拉法得到的系统状态方程，结合后向欧拉法（其中h 为步长）•+++=11m m m hX X X 即可得到系统的差分方程3梯形法由前面的系统状态方程，结合梯形法)(211+••+++=m m m m X X h X X 即可得到系统的差分方程4二阶显式Adams 方法由前面的状态方程，结合二阶显式Adams 方法)51623(12211--++-+=m m m m m F F F h X X 即可得到系统的差分方程但是二阶显式Adams 法不能自起步，要使方程起步，需要知道开始的三个值，但是我们只知道第一个值。
经过分析后，二阶显式Adams 方法精度是二阶的，而梯形法精度也是二阶的，因此我们可以先借助梯形法得到输出的前三个值，以达到起步的目的，然后借助上面得到的差分方程对其进行求解。
二、实验波形：下图为前向欧拉法、后向欧拉法、梯形法、二阶显式Adams 方法的系统差分方程得到相应的输出波形：图1 h=410 时四种方法的输出波形图2 h=56-⨯时四种方法的输出波形10图3 h=510-时四种方法的输出波形图4 h=610-时四种方法的输出波形三、实验分析：由输出波形可以看到各种方法的特点（在图中蓝色线均表示连续系统模型的实际输出波形，红色线表示在对应方法下系统的输出波形。
）：1前向欧拉法和二阶显式Adams方法对步长的要求很强。
步长太大，最后的到的结果不是绝对收敛，而是发散。
在小步长下才显得收敛，这也从另一方面验证，步长越小，截断误差越小；2步长不能太小，太小的步长相应的舍入误差和累积误差也会增大；3前向欧拉法也可称为显式欧拉法，后向欧拉法也可称为隐式欧拉法，可以看到，后向欧拉法的稳定域要比前向欧拉法大，计算精度也要高一些。
电路计算机仿真实验报告
电路计算机仿真分析实验报告实验一直流电路工作点分析和直流扫描分析一、实验目的1、学习使用Pspice软件，熟悉它的工作流程，即绘制电路图、元件类别的选择及其参数的赋值、分析类型的建立及其参数的设置、Probe窗口的设置和分析的运行过程等。
2、学习使用Pspice进行直流工作点分析和直流扫描分析的操作步骤。
二、原理与说明对于电阻电路，可以用直观法（支路电流法、节点电压法、回路电流法）列写电路方程，求解电路中各个电压和电流。
PSPICE软件是采用节点电压法对电路进行分析的。
使用PSPICE软件进行电路的计算机辅助分析时，首先在capture环境下编辑电路，用PSPICE的元件符号库绘制电路图并进行编辑、存盘。
然后调用分析模块、选择分析类型，就可以“自动”进行电路分析了。
需要强调的是，PSPICE软件是采用节点电压法“自动”列写节点电压方程的，因此，在绘制电路图时，一定要有参考节点（即接地点）。
此外，一个元件为一条“支路”（branch）,要注意支路（也就是元件）的参考方向。
对于二端元件的参考方向定义为正端子指向负端子。
三、示例实验应用PSPICE求解图1-1所示电路个节点电压和各支路电流。
图1-1 直流电路分析电路图R2图1-2 仿真结果四、选做实验1、实验电路图（1）直流工作点分析，即求各节点电压和各元件电压和电流。
（2）直流扫描分析，即当电压源Us1的电压在0-12V之间变化时，求负载电阻R L中电流I RL随电压源Us1的变化曲线。
IPRINT图1-3 选做实验电路图2、仿真结果Is21Adc1.000AVs35Vdc3.200A R431.200A23.20VVs47Vdc1.200A 0VR142.800AIs32Adc 2.000A12Vdc2.800AIIPRINT3.200A10.60V 12.00V Is11Adc 1.000A18.80V 28.80V15.60V3.600VR222.800ARL13.200A18.80VVs210Vdc2.800A Is53Adc3.000AI42Adc图1-4 选做实验仿真结果3、直流扫描分析的输出波形图1-5 选做实验直流扫描分析的输出波形4、数据输出V_Vs1 I(V_PRINT2)0.000E+00 1.400E+00 1.000E+00 1.500E+00 2.000E+00 1.600E+00 3.000E+00 1.700E+00 4.000E+00 1.800E+00 5.000E+00 1.900E+00 6.000E+00 2.000E+00 7.000E+00 2.100E+00 8.000E+00 2.200E+009.000E+00 2.300E+001.000E+012.400E+001.100E+012.500E+001.200E+012.600E+00从图1-3可以得到IRL与USI的函数关系为：I RL=1.4+(1.2/12)U S1=1.4+0.1U S1 (公式1-1)五、思考题与讨论：1、根据图1-1、1-3及所得仿真结果验证基尔霍夫定律。
计算机仿真技术实验教案
课程教案课程名称：计算机仿真技术实验任课教师：汤群芳所属院部：电气与信息工程学院教学班级：电气1403-04班教学时间：2015—2016学年第2学期湖南工学院课程基本信息1 实验一熟悉MATLAB环境及基本运算（验证性实验）一、本次课主要内容1、熟悉MATLAB环境；2、掌握MATLAB常用命令；3、MATLAB变量与运算符。
二、实验目的与要求1、熟悉MATLAB开发环境；2、掌握矩阵、变量、表达式的各种基本运算。
三、教学重点难点重点：矩阵的运算；难点：无。
四、教学方法和手段课堂讲授、演示；巡回指导。
五、作业与习题布置完成实验报告2 实验一熟悉MATLAB环境及基本运算（验证性实验）一、实验目的1．熟悉MATLAB开发环境2．掌握矩阵、变量、表达式的各种基本运算二、实验原理1.熟悉MATLAB环境熟悉MATLAB桌面和命令窗口、命令历史窗口、帮助信息浏览器、工作空间浏览文件和搜索路径浏览器。
3.MATLAB变量与运算符变量命名规则如下：（1）变量名可以由英语字母、数字和下划线组成（2）变量名应以英文字母开头（3）长度不大于31个（4）区分大小写MATLAB中设置了一些特殊的变量与常量，列于下表。
表1 MATLAB的特殊变量与常量MATLAB运算符，通过下面几个表来说明MATLAB的各种常用运算符3表2 MATLAB算术运算符表3 MATLAB关系运算符表4 MATLAB逻辑运算符表5 MATLAB特殊运算4. MATLAB的一维、二维数组的访问45. MATLAB的基本运算表7 两种运算指令形式和实质内涵的异同表6.MATLAB的常用函数表8 标准数组生成函数5表9 数组操作函数7.多项式运算poly——产生特征多项式系数向量roots——求多项式的根p=poly2str(c,‘x’)—（将特征多项式系数向量c转换为以习惯方式显示是多项式）conv, convs——多项式乘运算deconv——多项式除运算polyder(p)——求p的微分polyder(a, b)——求多项式a,b乘积的微分[p,q]=polyder(p1,p2)——求解多项式p1/p2微分的有理分式Polyval(p,A)——按数组运算规则求多项式p在自变量A的值polyvalm(p,A)——按矩阵运算规则求多项式p在自变量A的值三、实验仪器设备与器材计算机（安装有MATLAB软件平台）。
计算机仿真试验指导书
计算机仿真技术实验指导书1计算机仿真技术实验指导书MATLAB是一种交互式的以矩阵为基本数据结构的系统。
在生成矩阵对象时，不要求明确的维数说明。
所谓交互式，是指MATLAB的草稿纸编程环境。
与C语言或FORTRON语言作科学数值计算的程序设计相比较，利用MATLAB可节省大量的编程时间。
本实验指导书主要讨论四个实验。
实验一信号与系统的时域分析以及信号合成与分解1. 实验目的(1) 连续时间信号的向量表示法和符号运算表示法，典型离散信号表示；(2) 连续信号和离散信号的时域运算与时域变换；(3) 连续系统和离散系统的卷积，以及冲激响应、阶跃响应、单位响应、零状态响应；(4) 周期信号的傅立叶级数分解与综合（以周期方波为例）；2. 实验原理与方法(1) 信号在MA TLAB中的表示方法MATLAB用两种方法来表示连续信号，一种是用向量的方法来表示信号，另一种则是符号运算的方法来表示信号。
用适当的MATLAB语句表示出信号后，就可以利用MATLAB的绘图命令绘制出直观的信号时域波形。
向量表示法表示信号的方法是：MATLAB用一个向量表示连续信号的时间范围，另一个向量表示连续信号在该时间范围内的对应样值。
如下列代码p=0.001;t=-pi:p:pi;f=1+cos(t);plot(t,f)title('f(t)=1+cos(t)')xlabel('t')axis([-pi,pi,-0.2,2.4])执行后即可绘制连续信号1+cos(t)的时域波形。
借助于符号运算以及符号绘图函数ezplot，也可以绘制连续信号时域波形。
如下列代码syms tf=sym('1+cos(t)') %定义符号表达式ezplot(f,[-pi,pi]) %绘制符号表达式波形set(gcf,'color','w') %设置当前图形背景颜色为白色执行后即可绘制连续信号1+cos(t)的时域波形。
计算机仿真实验实验报告4次
（3）模块复制及删除；
（4）模块的调整
（5）ቤተ መጻሕፍቲ ባይዱ块参数的设置
（6）模块的连接
3、系统仿真运行
（1）在Simulink模型窗口下仿真
（2）在MATLAB命令窗口下仿真
4、仿真结果的输出和保存
（1）利用Scope模块；
（2）利用Out模块（在sinks库中），数据保存在MATLAB工作空间中(有tout和yout两项，分别为仿真时间向量和仿真输出向量），供以后调用和分析；
系统零极点增益模型
状态空间模型
系统模型的转换
系统模型参数的获取
时间延迟系统建模
模型属性设置和获取
系统模型的连接。
实验题目：
1、见教材《机电系统动态仿真》（机械工业出版社，刘白雁编）P-104，习题3、4。
获取已建立模型的参数；
2、系统开环传递函数为：
绘制当K=5、30时系统的Bode图，并判断系统的稳定性；计算K=5、30时系统的幅值和相位裕度；绘制K=5、30时系统的Nyquist图；绘制K=5、30时系统的Nichols图。
转速调节器是调速系统的主导调节器，它使转速n很快地跟随给定电压变化，稳态时可减小转速误差，如果采用PI调节器，则可实现无静差。电流调节器是内环调节器，在外环转速的调节过程中，它的作用是使电流紧紧跟随外环调节器的输出量变化。
实验要求：
根据直流电动机双闭环调速系统机构图，可以建立系统的Simulink仿真模型，经过仿真后，对所得结论进行分析比较，提高系统的动态性能。
3.在[–6，2]范围内用plot和fplot函数分别绘制二维曲线图。
4．绘制z=sin(x)*cos(y)的三维网格和三维曲面图，x,y变化范围均为[0，2π]。
计算机仿真实验报告
计算机仿真实验报告《计算机仿真》上机实验报告姓名：学号：专业：班级：实验⼀常微分⽅程的求解及系统数学模型的转换⼀．实验⽬的通过实验熟悉计算机仿真中常⽤到的Matlab指令的使⽤⽅法，掌握常微分⽅程求解指令和模型表⽰及转换指令，为进⼀步从事有关仿真设计和研究⼯作打下基础。
⼆. 实验设备个⼈计算机，Matlab软件。
三. 实验准备预习本实验有关内容（如教材第2、3、5章中的相应指令说明和例题），编写本次仿真练习题的相应程序。
四. 实验内容1. Matlab中常微分⽅程求解指令的使⽤题⽬⼀：请⽤MATLAB的ODE45算法分别求解下列⼆个⽅程。
要求:1.编写出Matlab 仿真程序；2.画出⽅程解的图形并对图形进⾏简要分析；3.分析下列⼆个⽅程的关系。
1．2．1.仿真程序xp1=@(t,x)-x^2;[t,x]=ode45(xp1,[0,20],[1]);plot(t,x);title('xp1')gridxp2=@(t,x)x^2;[t,x]=ode45(xp2,[0,20],[-1]);figure (2);plot(t,x);xlabel('t');ylabel('x')grid2.⽅程解的图形图形进⾏简要分析3.⼆个⽅程的关系题⽬⼆：下⾯⽅程组⽤在⼈⼝动⼒学中，可以表达为单⼀化的捕⾷者-被捕⾷者模式（例如，狐狸和兔⼦）。
其中1x 表⽰被捕⾷者， 2x 表⽰捕⾷者。
如果被捕⾷者有⽆限的⾷物，并且不会出现捕⾷者。
于是有1'1x x ，则这个式⼦是以指数形式增长的。
⼤量的被捕⾷者将会使捕⾷者的数量增长；同样，越来越少的捕⾷者会使被捕⾷者的数量增长。
⽽且，⼈⼝数量也会增长。
请分别调⽤ODE45、ODE23算法求解下⾯⽅程组。
要求编写出Matlab 仿真程序、画出⽅程组解的图形并对图形进⾏分析和⽐较。
function xp3=xp3(t,x)xp3=[x(1)-0.1*x(1)*x(2)+0.01*t ;-x(2)+0.02*x(1)*x(2)+0.04*t] [t,x]=ode45('xp3',[0,20],[30;20]); plot(t,x);grid[t,x]=ode23('xp3',[0,20],[30;20]); plot(t,x);gridODE45ODE232. Matlab 中模型表⽰及模型转换指令的使⽤题⽬三：若给定系统的的传递函数为1132106126)(23423+++++++=s s s s s s s s G 请⽤MATLAB 编程求解其系统的极零点模型。
计算机仿真实验
计算机仿真实验计算机仿真实验各位同学：（1）这是计算机仿真实验内容，给出的实验过程已在matlab 上运行，你们可以直接copy 在matlab 上运行，观察运行结果，对比学习。
（2）实验要求写出实验报告，实验报告的内容就按照以下给出的内容来写。
实验一：MATLAB 基础入门实验目的：熟悉MATLAB 环境，掌握一维数组的创建，二维数组的创建。
(1)一维数组的创建：观察下列计算结果，理解数组运算的意义：a=[1,2,pi,9,0]b=0:2:10c=linspace(1,2,10)a(2) b(5) c(6)(2)二维数组的创建：a=[1,2,3,0;9,22,1,1];观察下列计算结果，理解数组运算的意义：a(1,2) a(2,3) a(:) a(:,:) a(:,1) a(2,:)[1,2;3,4]+10[1,2;3,4]*[0.1,0.2;0.3,0.4][1,2;3,4].\[20,10;9,2]sin([1,2;3,4])[a,b]=find([1,2;3,4]>=[4,3;2,1]) [a,b]find([1,2;3,4]>=[4,3;2,1])(3) 观察下列计算结果，理解这些命令的意义Clear a ； Which ；clc 的作用实验目的：掌握MA TLAB 中基本的二维绘图plot(x,y)；plot(x,y,string)指出以下各个绘图命令的输出图形分别是什么，并上机验证t=0:pi/20:2*pi;y=sin(t);plot(t,y)t=[0 1]; x=[1 2]; y=[x;3 4]; z=[y;5 6];plot(t,x,'r') plot(t,y, 'b') plot(t,y') plot(t,z) plot(t,z') 线型 -实线 :虚线 -. 点划线 --间断线点标记 . 点 o 小圆圈 x 叉子符 +加号 *星号s 方格 d 菱形 ^朝上三角 v 朝下三 >朝右三角颜色 y 黄色 m 品红色 c 青色 r 红色 g 绿色b 蓝色w 白色上机指出以下各个绘图命令的输出图形分别是什么t=0:pi/20:2*pi; plot(t,sin(t) ,'r:>')hold on;plot(t,cos(t),'b-.h')title('sin(t),cos(t)的函数图形')xlabel('t=0:pi/20:2*pi;')ylabel('sin(t),cos(t)')legend('sin(t)','cos(t)')体验grid on /grid off;hold on/hold off;figure(2)的作用利用 plot 函数在一个坐标系下绘制以下函数的图形：y1=sin(x),y2=cos(x),y3=sin(2*x),x在0到2*pi区间y 1 用黑色间断线点标记为星号y 2 用红色实线点标记为小方格y 3 用蓝色虚线点标记为小圆圈实验目的：掌握matlab 编程的基本知识。
计算机仿真实验教案
计算机仿真实验教案实验一MATLAB基本操作（验证性实验）一、本次课主要内容1．熟悉MATLAB基本的操作界面。
2．掌握MATLAB中变量、数组、向量和矩阵等对象的生成和基本基本运算方法。
二、教学目的与要求1．熟悉MATLAB基本的操作界面。
2．掌握MATLAB中变量、数组、向量和矩阵等对象的生成和基本基本运算方法。
三、教学重点难点MATLAB中变量、数组、向量和矩阵等对象的生成和基本基本运算方法。
四、教学方法和手段课堂讲授、提问、讨论；使用多媒体教学方式。
五、作业与习题布置思考题实验一MATLAB基本操作（验证性实验）一、实验目的1．熟悉MATLAB基本的操作界面。
2．掌握MATLAB中变量、数组、向量和矩阵等对象的生成和基本运算方法。
3．掌握MATLAB中绘图的基本操作。
4．掌握MATLAB中的常用帮助命令使用方法二、实验原理1．常见数学函数表1-1 常见数学函数2、系统的在线帮助（1）help 命令：1）当不知系统有何帮助内容时，可直接输入help以寻求帮助: >> help（回车）2）当想了解某一主题的内容时，如输入：>> help syntax （了解Matlab 的语法规定）3）当想了解某一具体的函数或命令的帮助信息时，如输入：>> help sqrt （了解函数sqrt 的相关信息）（2）lookfor 命令现需要完成某一具体操作，不知有何命令或函数可以完成，如输入：>> lookfor line （查找与直线、线性问题有关的函数）3．常量与变量系统的变量命名规则：变量名区分字母大小写；变量名必须以字母打头，其后可以是任意字母，数字，或下划线的组合。
此外，系统内部预先定义了几个有特殊意义和用途的变量，见表1-2：表1-2 特殊的变量、常量（1）数值型向量（矩阵）的输入1）任何矩阵（向量），可以直接按行方式输入每个元素：同一行中的元素用逗号（，）或者用空格符来分隔；行与行之间用分号（；）分隔。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《仿真技术与应用》实验报告计算机仿真技术实验报告实验三利用数值积分算法的仿真实验实验三 利用数值积分算法的仿真实验实验目的1) 熟悉MATLAB 勺工作环境；2) 掌握MATLAB 勺.M 文件编写规则，并在命令窗口调试和运行程序;3)掌握利用欧拉法、梯形法、二阶显式Adams 法及四阶龙格库塔法构建系统仿 真模型的方法，并对仿真结果进行分析。
实验内容系统电路如图2.1所示。
电路元件参数：直流电压源E =1V ，电阻R=10门，电感L =0.01H ，电容C 二WF 。
电路元件初始值：电感电流i L (0) =0A ,电容电压u c (0) = 0V 。
系统输出量为电容电压u c (t)。
连续系统输出响应u c (t)的解析解为:u c (t)二 U s (1 — e_at (cos ■ t si nt a/ J)三、要求1) 利用欧拉法、梯形法、二阶显式 Adams 法及显式四阶Runge-Kutta 法构建系统仿真模 型，并求出离散系统的输出量响应曲线；(2-1)其中, a =2L图2.1 RLC 串联电路2) 对比分析利用欧拉法、梯形法、二阶显式 Adams 法及显式四阶Runge-Kutta 法构建系 统仿真模型的仿真精度与模型运行的稳定性问题；3) 分别编写欧拉法、梯形法、二阶显式 Adams 法及显式四阶Runge-Kutta 法的.m 函数文 件，并存入磁盘中。
.m 函数文件要求输入参数为系统状态方程的系数矩阵、仿真时间及仿 真步长。
编写.m 命令文件，在该命令文件中调用已经编写完成的上述 .m 函数文件，完成 仿真实验；4) subplot 和plot 函数将输出结果画在同一个窗口中，每个子图加上对应的标题。
四. 实验原理(1) 连续系统解析解连续系统输出响应u c (t)的解析解为：u c (t)二U s (1-e^t (cos t si nt x/ ，))(2) 原系统的传递函数根据所示电路图，我们利用电路原理建立系统的传递函数模型，根据系统的传递函数 是在零初始条件下输出量的拉普拉斯变换与输入量的拉普拉斯变换之比，可得该系统的传 递函数：(3) 系统的仿真模型在连续系统的数字仿真算法中，较常用的有欧拉法、梯形法、二阶显式Adams 法及显式 四阶Runge-Kutta 法等。
欧拉法、梯形法和二阶显式 Adams 法是利用离散相似原理构造的 仿真算法，而显式四阶 Runge-Kutta 法是利用Taylor 级数匹配原理构造的仿真算法。
对于 线性系统，其状态方程表达式为：；X(t) = Ax (t) + Bu (t)W)二 Cx (t)Du (t)其中: x h X 't)X 2(t)…x n (t) T 是系统的n 维状态向量u (t) = U i (t)比⑴ …U m (t )T 是系统的m 维输入向量其中, G(s)二U c (s) E(s)1/ LC s 2 R/ Ls 1/LCX (t o) = X 0X 页，y(t)二y i(t) y2(t)…y r(t)T是系统的r维输出向量A为n n阶参数矩阵，又称动态矩阵，B为n m阶输入矩阵，C为r n阶输出矩阵, D为r m阶交联矩阵。
根据图所示电路，系统状态方程模型：x(t)二Ax(t) B E y(t)二ex(t)式中，状态变量x = [x1,x2]T二[i L,u C]T，输出变量y(t)=u C，系数矩阵为:A「― R/L -1/U D刁/L l 厂r A=| ，B=| ，C = 101」。
]1/C 0」]0 一(1) 欧拉法利用前向欧拉法构建线性系统的仿真模型为：x m 1 二x m * x m h 二—A h x m - h Bu m,y m 1 二Cx m -1 ' DU m -1式中，h为积分步长，1为单位矩阵。
利用后向欧拉法构建线性系统的仿真模型为：. -1X m ・1 二x m - xm』二1 - A h X g h BU m 1.y m ・1 二Cx m -1 ' Du m -1对于前向欧拉法，系数矩阵为：A z = 1 h A，B z二h B，C z二C，D=0。
对于后向欧拉法，系数矩阵为：£ = 1 -h A ' , B z二1 - h A 一1h B，C z二C。
(2) 梯形法利用梯形法构建线性系统的仿真模型为：对图所示的系统，利用梯形法构造的系统差分方程具有形式:X m 1二 A t X m 2B t Ey m二 C t X m -1二1其系数矩阵为：A = H - h A / 21 h A /2 ,,B t 二 1 - h A /2 J h B 2，C t 二 C ，D = 0。
(3)二阶显式Adams 法利用二阶显式Adams 法构建线性系统的仿真模型为：h --]I x 羽=x + 一 23F —16F q + 5F 21 m 12 - m m _1m -2』m+1 = C x m 41 + DU m+1F m = Ax m - BU m式中：Fm」二 A xmdBUmJFm 上二 Ax m/ BU m/二阶显式Adams 法为多步计算方法，利用多步计算方法对系统进行仿真时， 需要与之 具有相同计算精度的单步计算方法辅助计算。
二阶显式Adams 法的计算精度为二阶，可以 采用梯形法或改进的Euler 法等辅助计算。
利用改进的Euler 法构建线性系统的仿真模型为：k 1 = f (t m ，X m ) = Ax m + B E k 2 二 f(t mh ,x mhkj = A x m k [h B E其中，m = 0,1。
由式计算出 片和x 2后，便可以转入由二阶显式 Adams 法构造 的离散系统模型计算，即系统差分方程。
其计算方程为：xm 单 x m 1y^i匚h /r h jh 口 11 — 一 A < 2丿1 + — A < 2丿x m + ㊁ B (Um + U ^1 )号 k 1 k 2)x m号 X m - Xm .1C X m 1D U m 1x mF m = Ax mB EFm 二 Ax m _1 B EFm^= Ax mJ2'B Ey m 1 二 Cx m 1(4)显式四阶 Runge-Kutta 法利用显式四阶Runge-Kutta 法构建线性系统的仿真模型为:k [=f(t m,xm)=Axm+ BEk 2 = f (t mh + —2xm+ h k 1)=A ( X m + k-i h / 2) + B E k 3= f(t m+h, 2xm+h k 2)=A (X m+ k 2h / 2) + B Ek 4 = f (t m+ h, X m + hk3)= A (X m + k3h ) + B EX m41 =x m h + - 6(k 1 + 2k 2 + 2k 3 + k 4)[_y m 舟=CxmH!五. 实验过程1■实验程序(1) 前向欧拉法fun ctio n []=RLC(R 丄,C,U,t,h) R=10; L=0.01; C=1.0e-6; U=1; t=0.01; h = 2.0e-4; m = fix(t/h); n = 2;A = [-R/L -1/L;1/C 0];B = [1/L;0]; D = [0 1]; E = [1 0;0 1]; %前向欧拉法%X m -1 X m—23F m -16F125F m _2for i=1:1: nx1(1: n,1) = 0;endfor k=1:mx1(1: n,k+1) = x1(1: n,k) + (A* x1(1: n,k)+B)*h;endfor k=1:1:my1(k) = D*x1(1: n,k);end%解析解%p = R/(2*L);w=sqrt(1/(L*C)-(R/(2*L)F2);for k=1:1:my(k) = U*(1-exp(-p*(k-1)*h) * ( cos(w*(k-1)*h) + si n(w*(k-1)*h)*p/w)); end %输出曲线%for k=1:1:mt(k) = (k-1)*h;endsubplot(2,3,1),plot(t,y,'g',t,y1,'r')legend('y解析解,','y1前向欧拉')title('前向欧拉法')(2)后向欧拉法fun ctio n []=RLC(R ,L ,C,U,t,h)R=10;L=0.01;C=1.0e-6;U=1;t=0.01;h = 2.0e-4;m = fix(t/h);n = 2;A = [-R/L -1/L;1/C 0];B = [1/L;0];D = [0 1];E = [1 0;0 1];%后向欧拉法%for i=1:1: nx2(1: n,1) = 0;endA1 = in v(E-A*h);for k=1:mx2(1: n,k+1) = A1*(x2(1: n,k) + B*h);endy2(k) = D*x2(1: n,k);end%解析解%p = R/(2*L);w=sqrt(1/(L*C)-(R/(2*L)F2);for k=1:1:my(k) = U*(1-exp(-p*(k-1)*h) * ( cos(w*(k-1)*h) + sin(w*(k-1)*h)*p/w));end%输出曲线%for k=1:1:mt(k) = (k-1)*h;endsubplot(2,3,2),plot(t,y,'g',t,y2,'r')legend('y解析解,','y2后向欧拉') title('后向欧拉法')(3)梯形法fun ctio n []=RLC(R ,L ,C,U,t,h)R=10;L=0.01;C=1.0e-6;U=1;t=0.01;h = 2.0e-4;m = fix(t/h);n = 2;A = [-R/L -1/L;1/C 0];B = [1/L;0];D = [0 1];E = [1 0;0 1];%梯形法%for i=1:1: nx3(1: n,1) = 0;endA2 = in v(E-A*h/2);for k=1:mx3(1: n,k+1) = A2*( x3(1: n,k) + B*h + A*x3(1: n,k)*h/2); end for k=1:1:my3(k) = D*x3(1: n,k);end%解析解%p = R/(2*L);w=sqrt(1/(L*C)-(R/(2*L)F2);y(k) = U*(1-exp(-p*(k-1)*h) * ( cos(w*(k-1)*h) + sin(w*(k-1)*h)*p/w));end%输出曲线%for k=1:1:mt(k) = (k-1)*h;endsubplot(2,3,3),plot(t,y,'g',t,y3,'r')legend('y解析解,','y3梯形法')title('梯形法')(4)二阶显式Adams法fun ctio n []=RLC(R ,L ,C,U,t,h)R=10;L=0.01;C=1.0e-6;U=1;t=0.01;h = 2.0e-4;m = fix(t/h);n = 2;A = [-R/L -1/L;1/C 0];B = [1/L;0];D = [0 1];E = [1 0;0 1];% 二阶显示Adams法%for i=1:1: nx4(1: n,1) = 0;endfor k=1:mx4(1: n, k+1) = A2*(x4(1: n,k) + B*h + A*x4(1: n,k)*h/2);endfor k=3:mfm1 = 23*(A*x4(1: n,k)+ B);fm2 = -16*(A*x4(1: n,k-1)+ B);fm3 = 5*(A*x4(1: n,k-2)+ B);x4(1: n, k+1) = x4(1: n,k)+(fm1+fm2+fm3)*h/12;endfor k=1:1:my4(k) = D*x4(1: n,k);end%解析解%p = R/(2*L);w=sqrt(1/(L*C)-(R/(2*L)F2);y(k) = U*(1-exp(-p*(k-1)*h) * ( cos(w*(k-1)*h) + sin(w*(k-1)*h)*p/w));end%输出曲线%for k=1:1:mt(k) = (k-1)*h;endsubplot(2,3,4),plot(t,y,'g',t,y4,'r')legend('y解析解，','y4Adams 法') title('二阶显式Adams 法')(5)四阶Runge-Kutta 法fun ctio n []=RLC(R ,L ,C,U,t,h)R=10;L=0.01;C=1.0e-6;U=1;t=0.01;h = 2.0e-4;m = fix(t/h);n = 2;A = [-R/L -1/L;1/C 0];B = [1/L;0];D = [0 1];E = [1 0;0 1];% 四阶Runge-Kutta 法%for i=1:1: n % 状态变量初值x5(1: n,1) = 0;endfor k=1:mx5(1: n, k+1) = A2*( x5(1: n,k) + B*h + A*x5(1: n,k)*h/2);endfor k=1:1:mk1= A*x5(1: n,k+1);k2=A*(x5(1: n,k+1)+h*k1/2);k3=A*(x5(1: n,k+1)+h*k2/2);k4=A*(x5(1: n,k+1)+h*k3);x5(1: n,k+1)=x5(1: n,k+1)+h*(k1+2*k2+2*k3+k4)./6;endfor k=1:1:my5(k) = D*x5(1: n,k);end%解析解%《仿真技术与应用》实验报告-11 -p = R/(2*L); w=sqrt(1/(L*C)-(R/(2*L)F2); for k=1:1:my(k) = U*(1-exp(-p*(k-1)*h) * ( cos(w*(k-1)*h) + sin(w*(k-1)*h)*p/w)); end%输出曲线%for k=1:1:mt(k) = (k-1)*h;end subplot(2,3,5),plot(t,y,'g',t,y5,'r')legend('y 解析解,','y5Runge-Kutta 法') title('显式四阶 Runge-Kutta 法')2■仿真图形取积分步长h=2*10-4s ,可以得到以下几个仿真图形: (1)前向欧拉法1x 10-1-2-3-4-50 0.001 0.002 0.003i* 0.0040.0050.006r0.0070.0080.0090.01(2)后向欧拉法后向欧拉法1.6 1.4 1.2 1 0.8 0.6 0.4 0.2y 解析解，-y2后向欧拉“0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01y 解析解， y1前向欧拉《仿真技术与应用》实验报告-12-显式四阶Runge-Kutta 法1 X 10 16前向欧y%法析解y1前向欧拉后向欧拉法0 -1 22梯形法y 解析解，1.5y2后向欧拉x 10■■二 阶显式Adams 法y 解析解， -21y3梯形法1.5y 解析解， -3 1 0y4Adams法-40.5-11一-50.0050.01-20.5/1 _0 00.005-30.01一11I \0.005 0.01-41 \ F0 0.0020.0040.0060.008 0.01-0.5-10.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01(3)梯形法梯形法21.8 1.6 1.4 1.2 --------- y 解析解,--------- y 3梯形法1 0.8 0.6 0.4 0.2 00.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01(4)二阶显式Adams 法0.5 (26)X 10二 阶 显 式 Adamsy 解析解，y4Adams 法-0.5-1-1.5-2-2.5-3x 10161前向欧拉法y 解析解， y1前向欧拉后向欧拉法-1-2 -3-4-50.0050.011.50.50.0050.001 0.002 0.003 0.004 0.005 y 解析解， y2后向欧拉梯形法1.50.50.010.0050.006 0.007 0.008 y 解析解， y3梯形法0.010.009 0.01(5)四阶 Runge-Kutta 法 四阶 Runge-Kuttay 解析解，y5龙格库 塔------------- y 解析解，------------- y5龙格库塔-10 0.0051.5《仿真技术与应用》实验报告六.实验结论1. 从仿真的稳定性看，当选取不同的积分步长时，欧拉法稳定性最低，梯形法稳定性其次，而显式四阶Runge-Kutta法、二阶显示Adams法稳定性较好。