当前位置：文档之家 › 哈工大计算机仿真实验三利用数值积分算法的仿真实验

哈工大计算机仿真实验三利用数值积分算法的仿真实验

  • 格式：pdf
  • 大小：458.31 KB
  • 文档页数：10

实验三利用数值积分算法的仿真实验一、实验目的1) 熟悉MATLAB的工作环境；2) 掌握MATLAB的 .M文件编写规则，并在命令窗口调试和运行程序；3) 掌握利用欧拉法、梯形法、二阶显式Adams法及四阶龙格库塔法构建系统仿真模型的方法,并对仿真结果进行分析。

二、实验内容与要求系统电路如图3所示。

电路元件参数：直流电压源，电阻，电感，电容。

电路元件初始值：电感电流，电容电压。

系统输出量为电容电压。

试利用欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法构建系统仿真模型，并求出离散系统的输出量响应曲线。

连续系统输出响应的解析解为：其中，三、实验原理在连续系统的数字仿真算法中,较常用的有欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法等。

欧拉法、梯形法和二阶显式Adams法是利用离散相似原理构造的仿真算法，而显式四阶Runge-Kutta法是利用Taylor级数匹配原理构造的仿真算法。

对于线性系统，其状态方程表达式为:(3-1)式（3-1）中，是系统的n维状态向量,是系统的m维输入向量，是系统的r维输出向量。

A为阶参数矩阵，又称动态矩阵，B为阶输入矩阵，C为阶输出矩阵，D为阶交联矩阵。

利用前向欧拉法构建线性系统的仿真模型为：（3-2）式中，为积分步长，为单位矩阵。

利用后向欧拉法构建线性系统的仿真模型为：（3-3）利用梯形法构建线性系统的仿真模型为：（3-4）利用二阶显式Adams法构建线性系统的仿真模型为：（3-5）式中： （3-6）二阶显式Adams法为多步计算方法，利用多步计算方法对系统进行仿真时，需要与之具有相同计算精度的单步计算方法辅助计算。

二阶显式Adams法的计算精度为二阶，可以采用梯形法或改进的Euler法等辅助计算。

利用改进的Euler法构建线性系统的仿真模型为：（3-7）利用显式四阶Runge-Kutta法构建线性系统的仿真模型为：（3-8）四、实验步骤与方法1. 建立系统数学模型根据图3所示电路，系统状态方程模型:(3-9)式中，状态变量，输出变量，系数矩阵为：，，。

2. 建立系统离散数学模型根据系统状态方程模型，利用式（3-2）~（3-8）给出的欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法的公式，分别构造与连续系统相似的离散系统模型，即系统的差分方程。

该模型是编写计算机仿真程序的基础。

对于图3所示的系统，利用欧拉法构造的系统差分方程具有形式：（3-10）对于前向欧拉法，式（3-10）的系数矩阵为：，，。

对于后向欧拉法，式（3-10）的系数矩阵为：，，。

对于图3所示的系统，利用梯形法构造的系统差分方程具有形式：（3-11）其系数矩阵为：，，。

利用二阶显式Adams法构造与连续系统相似的离散系统模型时，首先选择起步计算方法。

这里选择改进的Euler法。

其离散系统模型为：（3-12）其中，。

由式（3-12）计算出和后，便可以转入由二阶显式Adams法构造的离散系统模型计算，即系统差分方程。

其计算方程为：（3-13）（3-14）其中，。

利用显式四阶Runge-Kutta法构建的图3所示线性系统的离散系统模型为：（3-8）3. 编写Matlab的M文件程序利用欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法建立的离散系统差分方程编写系统仿真程序。

Matlab的M函数编写及运行见附录。

编写的M文件程序应能够满足实验要求。

4. 仿真实验为了比较分析算法特性，如稳定性、精度及误差，与积分步长的关系，应选择一组合理的积分对系统进行仿真实验。

试采用表3-1数据，分别利用欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法对图3电路进行仿真，给出仿真试验曲线及算法误差曲线（与准确解相比）。

从算法的稳定性、精度和算法误差与与积分步长的关系角度，对算法的仿真结果进行对比分析。

表3-1 积分步长的选择仿真时间t=0.01s积分步长/s五、实验报告要求1）给出欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法构建离散系统模型的过程；2）给出系统仿真程序；3）给出仿真试验曲线；4）分析仿真结果，从仿真模型实现的难易性、模型的稳定性、模型的精度等方面，对欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法构造的离散系统模型进行对比分析，并给出分析结论。

六、系统仿真程序function []=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: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% 后向欧拉法 %x2(1:n,1) = 0;endA1 = inv(E-A*h);for k=1:mx2(1:n,k+1) = A1*(x2(1:n,k) + B*h);endfor k=1:1:my2(k) = D*x2(1:n,k);end% 梯形法 %for i=1:1:nx3(1:n,1) = 0;endA2 = inv(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);endfor k=1:1:my3(k) = D*x3(1:n,k);end% 二阶显示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;endy4(k) = D*x4(1:n,k);end% 四阶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% 解析解 %p = R/(2*L);w=sqrt(1/(L*C)-(R/(2*L))^2);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,1),plot(t,y,'r',t,y1,'b')legend('y解析解,','y1前向欧拉')title('前向欧拉法')subplot(2,3,2),plot(t,y,'r',t,y2,'b')legend('y解析解,','y2后向欧拉')title('后向欧拉法')subplot(2,3,3),plot(t,y,'r',t,y3,'b') legend('y解析解,','y3梯形法')title('梯形法')subplot(2,3,4),plot(t,y,'r',t,y4,'b') legend('y解析解,','y4二阶显式Adams法') title('二阶显式Adams法')subplot(2,3,5),plot(t,y,'r',t,y5,'b') legend('y解析解,','y5显式四阶Runge-Kutta 法')title('显式四阶Runge-Kutta法')六、仿真实验曲线七、实验结论1）欧拉法为单步计算法，为自动起步的。

梯形法也为单步法，可以自动起步计算。

二阶显示Adams法需要知道k个初始值，不能自起步，二次函数很复杂，常用单步法来得到所需的几个初始值，很为复杂。

显式四阶Runge-Kutta法建模最为复杂，仿真时间也较长。

2）我们看到欧拉法最先与解析解分离开来，最后振荡开来，故稳定性最低；显式四阶Runge-Kutta法、二阶显示Adams 法和梯形法而则在较长时间可与曲线拟合收敛，故稳定性好。

3）显式四阶Runge-Kutta法，二阶显示Adams法模型的精度较高，函数曲线与真实曲线较为接近。

其次精确的是梯形法，取梯形面积，误差也较小。

前向欧拉法和后向欧拉法模型的精度最低。

4）二阶显示Adams法和梯形法对离散时间间隔要求低，前向欧拉法和后向欧拉法由于取的是矩形面积，离散时间间隔要求高。

5）从以上几种方法的仿真，共同规律是在小步长下都收敛，步长越小，误差越小，函数曲线与真实曲线较为接近；步长太大，最后的到的结果不是绝对收敛。

同时，同一算法隐式算法要比显式算法稳定性好，计算精度高。

合集下载

计算机仿真实验报告-实验三

计算机仿真实验报告-实验三

一、实验内容：实验三 利用欧拉法、梯形数法和二阶显式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、实验目的利用数学软件绘制数学函数曲线及曲面，通过实验了解函数图形的绘制方法。

2、实验内容⑴在同一个图形中，绘制双曲线，以及的双曲线2条渐近线。

⑵在同一个图形中，绘制球面与锥面相交的曲面。

⑶自选题目：绘制一个或者多个平面图形、空间曲面图形。

3、程序设计及运行结果(1)>> x=-5:0.1:5;ezplot('x^2-y^2=1');y1=x;y2=-x;hold on;plot(xy1);hold on;plot(xy2);(2) >> x=-5:0.1:5;y=x;z=x;[xyz]=meshgrid(xyz);f1=x.^2+y.^2+z.^2-1;f2=x.^2+y.^2-z;p1=patch(isosurface(xyzf10));set(p1 'FaceColor' 'm');p2=patch(isosurface(xyzf20));set(p2 'FaceColor' 'w');(3)>> x=-5:0.1:5;y=x;z=x;[xyz]=meshgrid(xyz);f1=x.^2+y.^2+z.^2-9;f2=x.^2+y.^2-2*z;p1=patch(isosurface(xyzf10));set(p1 'FaceColor' 'm');p2=patch(isosurface(xyzf20));set(p2 'FaceColor' 'w');4、讨论与分析在本次试验中初步了解了matlab，学会了一些简单绘图，加深了对函数的理解为以后实验作个铺垫，由浅入深的了解matlab.实验二：微积分实验1、实验目的熟悉并了解使用数学软件，进行微积分问题计算的相关数学软件命令，让学生通过实验理解微积分，解决微积分计算上的问题。

数值积分算法实验报告

数值积分算法实验报告

数值积分算法实验报告实验目的：验证数值分析理论实验工具：matlab摘要：本实验通过对比龙贝格积分算法和三点，五点勒让德高斯求积公式，对数值分析理论进行验证。

文章中提供了用matlab写的小程序即运行结果。

我们可以看到，龙贝格积分算法方便快捷且效率较高。

三点积分算法不足以满足精度要求，五点法基本可以满足，但是如果追求更高精度，则五点法改进困难，而龙贝格则可适用于任何精度的计算。

高斯求积的效率是比较高的，只需五点就可以达到非常高的精度。

一matlab源程序1．%龙贝格积分算法，jd表示相对精度function z=lbgjf(a,b,jd)h=(b-a);TT(1,1)=h.*(f(b)+f(a))/2;k=2;TT(1,2)=TT(1,1)./2+h/2.*f(a+h/2);TT(2,1)=TT(1,2).*4/3-TT(1,1)./3;z=TT(2,1);while abs((TT(k,1)-TT(k-1,1))./TT(k,1))>=jdk=k+1;h=h./2;for j=1:2.^(k-2)ff(1,j)=f(a+h*(j-1/2));endfff=sum(ff).*h/2;TT(1,k)=TT(1,k-1)./2+fff;for j=2:kTT(j,k-j+1)=4^(j-1).*TT(j-1,k-j+2)./(4^(j-1)-1)-TT(j-1,k-j+1)/(4^(j-1 )-1);z=TT(j,k-j+1);endend2 %五点法·¨function z=fivedlrd(a,b)Ak=[0.2369269 0.4786287 0.5688889 0.4786287 0.2369269];xk=[-0.9061798 -0.5384693 0 0.5384693 0.9061798];for i=1:5ff(i)=Ak(i).*f((b-a).*xk(i)./2+(b+a)./2);endz=(b-a)./2.*sum(ff)3%三点法¨function z=threedlrd(a,b)Ak=[0.5555556 0.8888889 0.5555556];xk=[-0.7745967 0 0.7745967];for i=1:3ff(i)=Ak(i).*f((b-a).*xk(i)./2+(b+a)./2); endz=(b-a)./2.*sum(ff)4%另用一matlab文件来表示函数function y=f(x)y=sin(x).^10;end二实验过程1 f(x)=x.^2,2 f(x)=x.^3可以看到，对于这种简单函数三种方法基本没有区别3 f(x)=sin(x).^10.lbgjf(2,3,10e-6)ans =0.0608>> fivedlrd(2,3)z =0.0608ans =0.0608>> threedlrd(2,3)z =0.0618ans =0.06184f(x)= exp(exp(sin(x).^10))threedlrd(2,3)z =2.9278ans =2.9278>> fivedlrd(2,3)z =2.9304ans =2.9304>> lbgjf(2,3,10e-6)ans =2.9304>>可见三点法已经不再满足要求。

哈工大 计算机仿真技术实验报告 仿真实验四基于Simulink控制系统仿真与综合设计

哈工大 计算机仿真技术实验报告 仿真实验四基于Simulink控制系统仿真与综合设计

基于Simulink 控制系统仿真与综合设计一、实验目的(1) 熟悉Simulink 的工作环境及其功能模块库； (2) 掌握Simulink 的系统建模和仿真方法；(3) 掌握Simulink 仿真数据的输出方法与数据处理；(4) 掌握利用Simulink 进行控制系统的时域仿真分析与综合设计方法； (5) 掌握利用 Simulink 对控制系统的时域与频域性能指标分析方法。

二、实验内容图2.1为单位负反馈系统。

分别求出当输入信号为阶跃函数信号)(1)(t t r =、斜坡函数信号t t r =)(和抛物线函数信号2/)(2t t r =时，系统输出响应)(t y 及误差信号)(t e 曲线。

若要求系统动态性能指标满足如下条件：a) 动态过程响应时间s t s 5.2≤；b) 动态过程响应上升时间s t p 1≤；c) 系统最大超调量%10≤p σ。

按图1.2所示系统设计PID 调节器参数。

图2.1 单位反馈控制系统框图图2.2 综合设计控制系统框图三、实验要求(1) 采用Simulink系统建模与系统仿真方法，完成仿真实验；(2) 利用Simulink中的Scope模块观察仿真结果，并从中分析系统时域性能指标（系统阶跃响应过渡过程时间，系统响应上升时间，系统响应振荡次数，系统最大超调量和系统稳态误差）；(3) 利用Simulink中Signal Constraint模块对图2.2系统的PID参数进行综合设计，以确定其参数；(4) 对系统综合设计前后的主要性能指标进行对比分析，并给出PID参数的改变对闭环系统性能指标的影响。

四、实验步骤与方法4.1时域仿真分析实验步骤与方法在Simulink仿真环境中，打开simulink库，找出相应的单元部件模型，并拖至打开的模型窗口中，构造自己需要的仿真模型。

根据图2.1 所示的单位反馈控制系统框图建立其仿真模型，并对各个单元部件模型的参数进行设定。

所做出的仿真电路图如图4.1.1所示。