matlab计算机仿真实验三.docx
《matlab计算机仿真实验三.docx》由会员分享，可在线阅读，更多相关《matlab计算机仿真实验三.docx（11页珍藏版）》请在冰点文库上搜索。
matlab计算机仿真实验三
实验3利用数值积分算法的仿真实验
一．实验目的
1）熟悉MATLAB的工作环境；
2）掌握MATLAB的.M文件编写规则，并在命令窗口调试和运行程序；
3）掌握利用欧拉法、梯形法、二阶显式Adams法及四阶龙格库塔法构建系统仿真模型的方法,并对仿真结果进行分析。
二．实验内容
系统电路如图2.1所示。
电路元件参数：
直流电压源
，电阻
，电感
，电容
。
电路元件初始值：
电感电流
，电容电压
。
系统输出量为电容电压
。
连续系统输出响应
的解析解为：
（2-1）
其中，
，
。
3．实验步骤
3.1求解各方法具体递推方程
3.1.1欧拉法
经分析得微分方程：
取状态变量：
则状态方程为：
由前向欧拉递推公式
得：
由后向欧拉递推公式
得：
解得：
3.1.2梯形法
由梯形法推导公式
得：
解得：
3.1.3二阶显式Adams法
由二阶显式Adams法递推公式
得：
而此算法为多步法，无法自动起步，所以需要补充同等精度的起步初始值。
此处采用梯形法求取。
3.1.4四阶龙格库塔法
由递推公式
得：
3.2用MATLAB仿真
1）根据离散方程和精确方程编写MATLAB程序,程序段如下：
E=1;
R=10;
L=0.01;
C=1e-6;
h=1e-6;
N=1e4;
x1
（1）=0;
x2
（1）=0;
a=R/（2*L）;
w=sqrt（1/（L*C）-（R/（2*L））^2）;
fork=1:
1:
N
x1（k+1）=x1（k）+h*x2（k）;
x2（k+1）=x2（k）+h*（1/（L*C））*（-R*C*x2（k）-x1（k）+E）;%前向欧拉法
uc（k）=E*（1-exp（-a*（h*k））*（cos（w*h*k）+sin（w*h*k）*a/w））;%解析解
end
t1=0:
h:
（h*N）;
t2=0:
h:
（h*（N-1））;
subplot（2,3,1）;
plot（t1,x1,'b',t2,uc,'r'）;
title（'前向欧拉法'）;
legend（'前向欧拉法','解析解'）;
fork=1:
1:
N
x2（k+1）=（L*C*x2（k）-h*x1（k）+h*E）/（h^2+h*R*C+L*C）;
x1（k+1）=x1（k）+h*x2（k+1）;%后向欧拉法
end
subplot（2,3,2）;
plot（t1,x1,'b',t2,uc,'r'）;
title（'后向欧拉法'）;
legend（'后向欧拉法','解析解'）;
fork=1:
1:
N
x1（k+1）=（（4*L*C+2*h*R*C-h^2）*x1（k）+4*h*L*C*x2（k）+2*h^2*E）/（h^2+2*h*R*C+4*L*C）;
x2（k+1）=2*（x1（k+1）-x1（k））/h-x2（k）;%梯形法
end
subplot（2,3,3）;
plot（t1,x1,'b',t2,uc,'r'）;
title（'梯形法'）;
legend（'梯形法','解析解'）;
fork=1:
1:
2
x1（k+1）=（（4*L*C+2*h*R*C-h^2）*x1（k）+4*h*L*C*x2（k）+2*h^2*E）/（h^2+2*h*R*C+4*L*C）;
x2（k+1）=2*（x1（k+1）-x1（k））/h-x2（k）;%梯形法取起步初始值
end
fork=3:
1:
N
x1（k+1）=x1（k）+h/12*（23*x2（k）-16*x2（k-1）+5*x2（k-2））;
x2（k+1）=x2（k）+h/12/（L*C）*（23*（-R*C*x2（k）-x1（k）+E）-16*（-R*C*x2（k-1）-x1（k-1）+E）+5*（-R*C*x2（k-2）-x1（k-2）+E））;
%显式Adams法
end
subplot（2,3,4）;
plot（t1,x1,'b',t2,uc,'r'）;
title（'显式ADams法'）;
legend（'显式Adams法','解析解'）;
fork=1:
1:
N
k1=x2（k）;
g1=1/L/C*（-R*C*x2（k）-x1（k）+E）;
k2=x2（k）+h/2*g1;
g2=1/L/C*（-R*C*（x2（k）+h/2*g1）-（x1（k）+h/2*k1）+E）;
k3=x2（k）+h/2*g2;
g3=1/L/C*（-R*C*（x2（k）+h/2*g2）-（x1（k）+h/2*k2）+E）;
k4=x2（k）+h*g3;
g4=1/L/C*（-R*C*（x2（k）+h*g3）-（x1（k）+h*k3）+E）;
x1（k+1）=x1（k）+h/6*（k1+2*k2+2*k3+k4）;
x2（k+1）=x2（k）+h/6*（g1+2*g2+2*g3+g4）;%四阶龙格库塔法
end
subplot（2,3,5）;
plot（t1,x1,'b',t2,uc,'r'）;
title（'四阶龙格库塔'）;
legend（'四阶龙格库塔法','解析解'）;
2）不同步长下的仿真图像
T=1e-7,N=1e5
T=1e-6,N=1e4
T=1e-5,N=1e3
T=5e-5,N=2e2
T=1e-4,N=4e1
T=1e-7,N=1e5局部放大图
4．实验结果分析
难易性：
欧拉法为单步计算法，给定初始值就可以直接迭代递推，方法较简单。
梯形法同样为单步计算法，对于状态方程需要进行简单的关系式转换，属于间接迭代递推。
显式Adams法为多步法，不能自起步，需要通过其他单步法计算需要的初始值，计算较复杂。
显式四阶龙格库塔法虽为单步法，但建模最复杂，计算量极大。
稳定性：
通过不同步长下的仿真曲线，可以判断出对于稳定性，前向欧拉法<后向欧拉法<显式Adams法<梯形法<四阶龙格库塔。
精确性：
通过小步长下的局部放大图可观察出，对于精确性，前向欧拉法˜后向欧拉法<显式Adams法˜四阶龙格库塔<梯形法。
升级会员 