内容发布更新时间 : 2024/12/25 3:27:03星期一 下面是文章的全部内容请认真阅读。
5. 设计M文件计算:
x=0:0.1:10
sum??(xi?2xi)2i?0100当sum>1000时停止运算，并显示求和结果及计算次数。
clear
x=0:0.1:10; s=1; sumx=0;
while(sumx<=1000)
sumx=sumx+x(s)*x(s)-2*x(s); s=s+1; end s sumx
二、实验内容：
M文件设计：仿真时间0～1秒，计算步长h=0.1，初始值y(0)=1，对连续系统
dy?x?y 采用以下方法建立仿真模型： dx① 前向Euler: ② 后向Euler: ③ 梯形法： ④ 改进Euler:
⑤ 经典RK4：
yn?1?yn?hfn
yn?1?yn?hfn?1hyn?1?yn??fn?fn?1?2?k1?f(xm,y(xm))?k?f(x?h,y(x)?hk)mm1?2h?y?y(x)?(k1?k2)?m?1m2?????????k1?f(xm,y(xm))k2?f(xm?hh,ym?k1)22k3?f(xm?hh,ym?k2)22k4?f(xm?h,ym?hk3)ym?1?ym?h?k1?2k2?2k3?k4?6yn?yn?1?h?23fn?1?16fn?2?5fn?3?12
⑥ Adams法： ⑦ Milne法： y?hm?1?ym?13?fm?1?4fm?fm?1?
并与精确解y(x)?2ex?1?x的结果进行比较。
%qianxiangEuler x1(1)=0;h=0.1;y1(1)=1; for s=1:1:10
f1(s)=x1(s)+y1(s); y1(s+1)=y1(s)+h*f1(s); x1(s+1)=x1(s)+h; end
%houxiangEuler y2(1)=1; x2(1)=0; h=0.1;
for s=1:1:10
y2(s+1)=(y2(s)+h*(x2(s)+h))/(1-h); x2(s+1)=x2(s)+h; end
%gaijinEuler x4(1)=0; y4(1)=1; h=0.1;
for s=1:1:10 k1=x4(s)+y4(s);
k2=x4(s)+h+y4(s)+h*k1; y4(s+1)=y4(s)+(k1+k2)*h/2; x4(s+1)=x4(s)+h;end
%梯形法 x3(1)=0; y3(1)=1;
h=0.1;for s=1:1:10
y3(s+1)=(y3(s)+(2*x3(s)+y3(s)+h)*h/2)/(1-h/2); x3(s+1)=x3(s)+h; end %RK4 x5(1)=0; y5(1)=1; h=0.1;
for s=1:1:10
k1=x5(s)+y5(s);
k2=x5(s)+h/2+y5(s)+k1*h/2; k3=x5(s)+h/2+y5(s)+k2*h/2;
k4=x5(s)+h+y5(s)+k3*h; y5(s+1)=y(s)+(k1+2*k2+2*k3+k4)*h/6; x5(s+1)=x5(s)+h;
2
end
-ams·¨ x6(1)=0; y6(1)=1; h=0.1;
y6(2)=y6(1)+h*(x6(1)+y6(1)); x6(2)=x6(1)+h;
y6(2)=y6(1)+(x6(1)+y6(1)+x6(2)+y6(2))*h/2; y6(3)=y6(2)+h*(x6(2)+y6(2)); x6(3)=x6(2)+h;
y6(3)=y6(2)+(x6(2)+y6(2)+x6(3)+y6(3))*h/2; for s=3:1:10
x6(s+1)=x6(s)+h;
y6(s+1)=y6(s)+(23*(y6(s)+x6(s))-16*(y6(s-1)+x6(s-1))+5*(y6(s-2)+x6(s-2)))*h/12; end
%Milne·¨ y7(1)=1; x7(1)=0; h=0.1;
y7(2)=y7(1)+h*(x7(1)+y7(1)); x7(2)=x7(1)+h;
y7(2)=y7(1)+(x7(1)+y7(1)+x7(2)+y7(2))*h/2; y7(3)=y7(2)+h*(x7(2)+y7(2)); x7(3)=x7(2)+h;
y7(3)=y7(2)+(x7(2)+y7(2)+x7(3)+y7(3))*h/2; for s=3:1:10
y7(s+1)=y7(s)+h*(x7(s)+y7(s)); x7(s+1)=x7(s)+h;
y7(s+1)=y7(s)+(x7(s)+y7(s)+x7(s+1)+y7(s+1))*h/2;
y7(s+1)=y7(s-1)+((y7(s+1)+x7(s+1))+4*(y7(s)+x7(s))+(y7(s-1)+x7(s-1)))*h/3; end
%jingquejie x=0:0.1:1;
y=2*exp(x)-1-x;
没有画图，画图部分自己加
1、函数文件设计：
设计一个函数文件实现一个阶乘运算 y(n)?n!?1?2?3???(n?1)?n， 为保证函数的通用性，当输入负数或小数时，显示出错提示： disp('Input parameter must be a positive integer!') 提示：fix(x) 对零方向取整数 ceil(x) 对+?方向取整数 round(x) 四舍五入取整数 并设计程序调用该函数。
function[y]=jc(n)%阶乘 if n<0
disp('Input parameter must be a positive integer!'); disp('亲 输入整数哦');
3