武汉理工大学过程仿真与建模作业

实实验验一一面面向向方方程程的的数数值值积积分分方方法法仿仿真真一、实验目的：加深理解四阶龙格---库塔法的原理及其稳定性二、实验内容对下列系统进行仿真A．线性定常系统'24(0)1u(01)tB．非线性系统()()()()()()()()dxtrxtaxtytdtdytsytbxtytdt其其中中：：1.r=0.001,a=2*104,s=0.015,b=10-4;x(0)=1200,y(0)=6002.r=0.001,a=2*10-6,s=0.01,b=10-6;x(0)=12000,y(0)=6000三三、、实实验验要要求求：：1.为保证稳定性，分析系统（1）的最大仿真步长（方法自选）2.设计MatLab、Fortran或C程序，用四阶龙格---库塔法进行仿真计算，改变参数及仿真步长，观察实验结果，寻找最合宜的仿真步长和临界仿真步长；四、实验报告1.实验所用程序清单2.实验结果及分析1.分析系统A的最大仿真步长，寻找临界仿真步长1.1四阶阶龙格—库塔的稳定性四阶龙格—库塔格式为：11234(22)6nnhuuKKKK其中1(,)nnKftu2112(,)2nnhKftuK3122(,)2nnhKftuK413(,)nnKftuhK由于'24uu，则有2233411(63)64nnnhuuhhhu为保证龙格—库塔的稳定性，则有1nnuu，从而可以得223341(63)64nnnhuhhhuu1.2通过编程求最大稳定步长clearallsymsK1K2K3K4yhy1=-24*y;y1=inline(y1);K1=y1(y);K2=y1(y+1/2*K1*h);K3=y1(y+1/2*K2*h);K4=y1(y+h*K3);R1=K1+2*K2+2*K3+K4;R11=inline(R1);R=solve(R1,'h');R=eval(R);p=R==real(R);R=R(p);fprintf('保证稳定性的最大仿真步长为%0.5f\n',R)结果显示为：1.3改变参数及仿真步长，观察实验结果MATLAB编写程序作图以步长h=0.025为例：clearallsymsK1K2K3K4yhy1=-24*y;y1=inline(y1);K1=y1(y);K2=y1(y+1/2*K1*h);K3=y1(y+1/2*K2*h);K4=y1(y+h*K3);R1=K1+2*K2+2*K3+K4;R11=inline(R1);h=0.025f=1;forii=1:10t=0.1*ii;tt(ii)=t;n=(t-0)/h;fori=1:nf=f+1/6*h*R11(h,f);endff(ii)=f;endplot(tt,ff,'o-')holdona=legend('步长h为0.025',1)xlabel('t')ylabel('u')1.31寻找临界步长分别取步长h为0.025、0.05、0.075、0.01、0.0125，用MATLAB作出图像0.10.20.30.40.50.60.70.80.9100.010.020.030.040.050.060.070.080.090.1tu步长h为0.025图1,步长h为0.025时，u的值0.10.20.30.40.50.60.70.80.9100.020.040.060.080.10.12tu步长h为0.05图2,步长h为0.05时，u的值0.10.20.30.40.50.60.70.80.9100.050.10.150.20.250.30.35tu步长h为0.075图3,步长h为0.075时，u的值0.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.7tu步长h为0.1图4,步长h为0.1时，u的值0.10.20.30.40.50.60.70.80.9100.511.522.533.5x105tu步长h为0.125图5,步长h为0.125时，u的值由图1-5可知，临界步长在0.1~0.125之间为了方便分析，计算A方程的真实值，并作图clearallt=0.1:0.1:1u=exp(-24*t)plot(t,u,'o-')legend('真实值',1)xlabel('t')ylabel('u')运行得到图6：0.10.20.30.40.50.60.70.80.9100.010.020.030.040.050.060.070.080.090.1tu真实值图6,u的真实值在0.1~0.125附近步长，得到图7如下：图7比较真实值和0.1附近步长可以得出结论，最合适仿真步长如果太大，结果不精确。如果太小，虽然结果精确，但是计算复杂。因此：最合适仿真步长取0.1。临界仿真步长为：0.11604。2.分析系统B的临界仿真步长和最合适步长2.1