计算方法上机实验试题

数学实验期末上机考核学号姓名年级专业学号姓名年级专业学号姓名年级专业学号姓名年级专业学号姓名年级专业请按要求回答以下问题：Part A、Matlab软件的应用（本大题两大类，共70分）以下问题的解答过程中如需软件，均应用matlab软件进行。

一、计算题（本大题共4小题，共55分）1、缆车的缆绳绳长问题（20分）某旅游景点从山脚到山顶有一缆车索道，全长约1471m，高差为380m，采用循环单线式修建，缆绳悬挂在下站到上站的行程的8个铁塔上，这8个铁塔依山势走向而距离不等（示意图见图1）。

假设从下站到第1铁塔的水平距离为d0，高差为h0，记此为第0段缆绳；从第1铁塔到第2铁塔的水平距离为d1，高差为h1，此为第1段缆绳；...，从第8铁塔到上站的水平距离为d8，高差为h8，此为第8段缆绳。

具体数据见下表：第0段第1段第2段第3段第4段第5段第6段第7段第8段d0 d1 d2 d3 d4 d5 d6 d7 d8 220 200 140 120 100 120 140 200 220 h0 h1 h2 h3 h4 h5 h6 h7 h850 45 40 38 34 38 40 45 50 假设每一段缆绳下垂的最低点不低于两端铁塔最低塔顶悬挂绳处1m，并假定所悬挂的缆绳为抛物线（图2）。

请按以下要求回答问题：图1 从山脚到山顶的缆车索道示意图图2 第i段所悬挂的抛物线缆绳（i=0,1, (8)（1）对第i段缆绳建立如下图所示的坐标系，请写出该段抛物线方程及所满足的条件。

（4分）（2）请计算第2段缆绳的方程各参数，及该段缆绳的长度，并画出该段缆绳的抛物线图形（10分）。

该问题的解答步骤如下：步骤一、编写程序计算第2段缆绳方程的各参数提示1：方程（组）求解命令：solve提示2：字符型与数值型的转换命令：double，例如double(2/5)=0.4步骤二、根据程序运行的结果，写出抛物线方程的数学表达式如下：步骤三、应用matlab 计算缆绳长度，请写明编程语言及计算结果。

数值计算方法I 上机实验考试题（两题任选一题）
１．小型火箭初始质量为900千克，其中包括600千克燃料。

火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。

设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）．重力加速度取9.8米/秒2.
A. 建立火箭升空过程的数学模型（微分方程）；
B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度．
２．小型火箭初始质量为1200千克，其中包括900千克燃料。

火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生40000牛顿的恒定推力．当燃料用尽时引擎关闭。

设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数记作k ，火箭升空过程的数学模型为
0)0(,0,01222==≤≤-+⎪⎭
⎫ ⎝⎛-==t dt dx x t t mg T dt dx k dt x d m 其中)(t x 为火箭在时刻t 的高度，m =1200-15t 为火箭在时刻t 的质量，T （=30000牛顿）为推力，g (=9.8米/秒2)为重力加速度, t 1 (=900/15=60秒)为引擎关闭时刻．
今测得一组数据如下（t ~时间（秒），x ~高度（米），v ~速度（米/秒））：
现有两种估计比例系数k 的方法：
1．用每一个数据(t,x,v )计算一个k 的估计值（共11个），再用它们来估计k 。

2．用这组数据拟合一个k ．
请你分别用这两种方法给出k 的估计值，对方法进行评价，并且回答，能否认为空气阻力系数k=0.5（说明理由）．。

（完整版）数值计算⽅法上机实习题答案1．设?+=105dx xx I nn ，（1）由递推公式nI I n n 151+-=-，从0I 的⼏个近似值出发，计算20I ；解：易得：0I =ln6-ln5=0.1823, 程序为：I=0.182; for n=1:20I=(-5)*I+1/n; end I输出结果为：20I = -3.0666e+010 （2）粗糙估计20I ，⽤nI I n n 515111+-=--，计算0I ；因为 0095.056 0079.01020201020≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2120=+=I 程序为：I=0.0087; for n=1:20I=(-1/5)*I+1/(5*n); end I0I = 0.0083（3）分析结果的可靠性及产⽣此现象的原因(重点分析原因)。

⾸先分析两种递推式的误差；设第⼀递推式中开始时的误差为000I I E '-=，递推过程的舍⼊误差不计。

并记nn n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。

因为=20E 20020)5(I E >>-，所此递推式不可靠。

⽽在第⼆种递推式中n n E E E )51(5110-==-=Λ，误差在缩⼩，所以此递推式是可靠的。

出现以上运⾏结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。

2．求⽅程0210=-+x e x的近似根，要求41105-+?<-k k x x ，并⽐较计算量。

（1）在[0，1]上⽤⼆分法；程序：a=0;b=1.0;while abs(b-a)>5*1e-4 c=(b+a)/2;if exp(c)+10*c-2>0 b=c; else a=c; end end c结果：c =0.0903（2）取初值00=x ，并⽤迭代1021x k e x -=+；程序：x=0; a=1;while abs(x-a)>5*1e-4 a=x;x=(2-exp(x))/10; end x结果：x =0.0905（3）加速迭代的结果；程序：x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;y=exp(x)+10*x-2; z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x); b=x; end x结果：x =0.0995（4）取初值00=x ，并⽤⽜顿迭代法；程序：x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x; end x结果： x =0.0905（5）分析绝对误差。

计算方法上机实习题目（四）——龙贝格算法计算椭圆周长一、 题目用龙贝格算法计算椭圆110040022=+y x 的周长，使误差不超过410-。

二、解题方法由题目中椭圆的标准方程 110040022=+y x 知，该椭圆的参数方程为)2,0[ sin 10cos 20π∈⎩⎨⎧==t ty t x 参考微积分中曲线弧长计算公式dt t y t x s ⎰+=βα)()(2'2' 知该椭圆的周长计算公式为 dt t t s ⎰+=π2022cos 100sin 400 故该问题为利用龙贝格算法计算数值积分。

参考教材式（5.42），若以一个二元数组T[i][j]代表式中的)(j i T ，则式（5.42）可化为⎪⎪⎪⎩⎪⎪⎪⎨⎧=⋯=--=--+-+=+-=-+-=-∑-)3,2,1;,2,1,0(144]2)12([221)]()([2][]1[]1[]1[][][21]1[]0[][]0[]0[]0[1m k T T T a b i a f a b T T b f a f a b T m k m k m m k m i l l l l l 其中t t t f 22cos 100sin 400)(+=，a ，b 为区间端点值，π2,0==b a 。

相应的计算程序段为T[0][0]=(b-a)*(f(a)+f(b))/2;k=1;T[0][1]=T[0][0]/2+(b-a)*Sum(1)/pow(2,1);T[1][0]=(4*T[0][1]-T[0][0])/(4-1);k=2;T[0][2]=T[0][1]/2+(b-a)*Sum(2)/pow(2,2);T[1][1]=(4*T[0][2]-T[0][1])/(4-1);T[2][0]=(pow(4,2)*T[1][1]-T[1][0])/(pow(4,2)-1);k=3;T[0][k]=T[0][k-1]/2+(b-a)*Sum(k)/pow(2,k);T[1][k-1]=(4*T[0][k]-T[0][k-1])/(4-1);T[2][k-2]=(pow(4,2)*T[1][k-1]-T[1][k-2])/(pow(4,2)-1);T[3][k-3]=(pow(4,3)*T[2][k-2]-T[2][k-3])/(pow(4,3)-1);for(k=4;fabs(T[3][k-4]-T[3][k-5])>=pow(10,-4);k++){T[0][k]=T[0][k-1]/2+(b-a)*Sum(k)/pow(2,k);T[1][k-1]=(4*T[0][k]-T[0][k-1])/(4-1);T[2][k-2]=(pow(4,2)*T[1][k-1]-T[1][k-2])/(pow(4,2)-1);T[3][k-3]=(pow(4,3)*T[2][k-2]-T[2][k-3])/(pow(4,3)-1);}其中f（x）和Sum（l）为调用的子函数，子函数程序如下：double Sum(int l){i;intdoublesum,a,b;double f(double x);a=0.0;b=2*PI;sum=0;for(i=1;i<=pow(2,l-1);i++)sum=sum+f(a+(2.0*i-1.0)*(b-a)/pow(2,l));sum;return}double f(double x){y;doubley=sqrt(pow(20*sin(x),2)+pow(10*cos(x),2));y;return}因为并不是一开始T[3][j]就同步出现，所以需将k=3之前的各值先单独算出，又因为要10-，所以控制循环结束的条件是fabs(T[3][k-4]-T[3][k-5])>=pow(10,-4)，求最后误差不超过4其中fabs（x）为求绝对值的函数，这里需要注意的是，每次判断时已经执行了k++这条指令，所以在判断是应该是T[3][k-4]-T[3][k-5]，而不是T[3][k-3]-T[3][k-4]，我在最初编写程序时就忽略了这个问题，导致花费很久调试程序。

《计算方法》上机实验试题
1. (25分)计算积分
dx x x I n n ⎰
+=105， n=0,1,2,…,20 若用下列两种算法 (A) n I I n n 151+
-=- (B) ⎪⎭
⎫ ⎝⎛-=-n n I n I 1511 试依据积分I n 的性质及数值结果说明何种算法更合理。

2. (25分)求解方程f(x)=0有如下牛顿迭代公式
)()(111---'-
=n n n n x f x f x x ， n ≥1，x 0给定 (1) 编制上述算法的通用程序，并以ε<--1n n x x （ε为预定精度）作为终止迭
代的准则。

(2)
利用上述算法求解方程 0cos :)(=-=x x x f
这里给定x 0=π/4，且预定精度ε=10-10。

3. (25分) 已知)3()(x x e x e x f -=，
(1)
利用插值节点x 0=1.00，x 1=1.02，x 2=1.04，x 3=1.06，构造三次Lagrange 插值公式，由此计算f(1.03)的近似值，并给出其实际误差； (2) 利用插值节点x 0=1，x 1=1.05构造三次Hermite 插值公式，由此计算f(1.03)的近
似值，并给出其实际误差。

4. (25分) 利用Romberg 算法计算积分
⎰
+4802cos 1dx x
精确到10
－4
总体要求：打印各题的程序代码及数值结果。

数值计算⽅法上机实习题考证--------------------------------------------------- 此⽂档包含我们计算⽅法的经典算法包含（数值计算⽅法上机实习题）1．设?+=105dx xx I nn ，（1）由递推公式n I I n n 151+-=-，从0I 的⼏个近似值出发，计算20I ；（2）粗糙估计20I ，⽤nI I n n 51511+-=-，计算0I ；（3）分析结果的可靠性及产⽣此现象的原因(重点分析原因)。

(1) 解答：n=0，0.1823)05ln()15ln()5(51515101010=+-+=++=+=+=x d xdx x dx x x I nn这⾥可以⽤for 循环，while 循环，根据个⼈喜好与习惯：for 循环程序： While 循环程序： I=0.1823; I=0.1823; for n=1:20 i=1;I=(-5)*I+1/n; while i<21 End I=(-5)*I+1/i; I i=i+1; fprintf('I20=%f',I) end I = -2.0558e+009 >> II20=-2055816073.851284>> I = -2.0558e+009 (2) 粗略估计I 20： Mathcad 计算结果： for 循环程序： While 循环程序： >> I=0.007998; I=0.007998; >> for n=1:20 n=1;I=(-0.2)*I+1/(5*n); while n<21End I=(-0.2)*I+1/(5*n); >> I n=n+1; I =0.0083 end >> II =0.0083(3) 算法误差分析：计算在递推过程中传递截断误差和舍⼊误差第⼀种算法：（从1——>20）1x x 205x +d 7.998103-?=*000e I I =-*115(5)5()555n n n n n n n n n n e I I I I I I e e e n n------=-=-+--+=-===误差放⼤了5n倍，算法稳定性很不好；第⼆种算法：（从20——>1）*n n ne I I =-***111111111()()555555n n n n n n nn e I I I I I I e n n ---=-=-+--+=-=0111...()55n ne e e ===误差在逐步缩⼩，算法趋近稳定，收敛。

数值计算方法上机实习题1． 设⎰+=105dx xx I nn ， （1） 由递推公式nI I n n 151+-=-，从0=0.1822I , 0=0.1823I 出发，计算20I ； （2） 20=0I ，20=10000I , 用nI I n n 51511+-=-，计算0I ；（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。

解：（1）程序如下: clear all clc I=0.1822; %题中的已知数据 for n=1:20; I=(-5)*I+1/n; %由递推公式所得 end fprintf('I20=%f\n',I) M=0.1823; %与I 的计算结果形成对比for i=1:20; M=(-5)*M+1/i; %由递推公式所得 end fprintf('M20=%f\n',M) 输出结果为： I20=-11592559237.912731 M20=-2055816073.851284 （2）程序如下: clear all clc I=0; %赋予I20的初始值 for n=0:19; I=(-1/5)*I+1/(5*(20-n)); %有递推公式得 end fprintf('I0=%f\n',I)M=10000; for i=0:19; M=(-1/5)*M+1/(5*(20-i));%有递推公式得 end fprintf('M0=%f\n',M) 输出结果为： I0=0.182322 M0=0.182322（3）由输出结果可看出第一种算法为不稳定算法，第二中算法为稳定算法。

由于误差*000***21111120115(5)5()555nn n n n n n n n n e I I e I I I I I I e e e n n------=-=-=-+--+=-===第一种算法为正向迭代算法，每计算一步误差增长5倍，虽然所给的初始值很接近，随着n 的增大，误差也越来越大。

目录1.计算方法A 上机作业 (1)上机练习目的 (1)上机练习任务 (1)计算方法A 上机题目 (1)程序设计要求 (1)上机报告要求 (1)2.QR 分解法求解线性方程组 (2)计算原理 (2)程序框图 (8)计算实习 (9)Matlab代码 (11)3.共轭梯度法求解线性方程组 (13)计算原理 (13)程序框图 (14)计算实习 (15)Matlab代码 (19)4.三次样条插值 (21)计算原理 (21)程序框图 (23)计算实习 (24)Matlab代码 (26)5.四阶龙格－库塔法求解常微分方程的初值问题 (30)计算原理 (30)程序框图 (31)计算实习 (32)Matlab代码 (33)1.计算方法A 上机作业上机练习目的❑ 复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。

❑ 利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。

上机练习任务•利用计算机语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。

•掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。

•写出上机练习报告。

计算方法A 上机题目1. QR 分解方法求解线性方程组。

（第二章）2. 共轭梯度法求解线性方程组。

（第三章）3. 三次样条插值（第四章）4. 四阶龙格－库塔法求解常微分方程的初值问题程序设计要求1. 程序要通用的，在程序设计时要充分考虑哪些变量应该可变的。

2. 程序要求调试通过。

上机报告要求报告容包括：● 每种方法的算法原理及程序框图。

● 程序使用说明。

● 算例计算结果。

2. QR 分解法求解线性方程组计算原理当n x R ∈是任意给定的非零向量，n v R ∈是任意给定的单位向量，则存在初等反射阵2T H I uu =-,使得Hx v σ=，其中σ为常数，当取单位向量2x vu x vσσ-=-时，由u 确定的矩阵H 必定满足Hx v σ=，所以在计算过程中取u 的值为上述值。

数值计算方法上机实验题数值计算方法上机实验实验内容：1.要求:分别用复化梯形法,复化Simpson 法和 Romberg 公式计算.2．给定积分dx e x31和dx x ?311 ,分别用下列方法计算积分值要求准确到510- ,并比较分析计算时间. 1）变步长梯形法; 2）变步长 Simpson 法; 3) Romberg 方法.算法描述：1、复合梯形法：?=tdt t a t V 0)()( ))()(2)((211∑-=++=n k k n b f x f a f hT输入被积函数数据点t,a. 输出积分值.n T复合Simpson 法：?=tdt t a t V 0)()( ))()(2)(4)((6101121∑∑---=++++=n k n k k k n b f x f x f a f h输入被积函数f(x),积分区间[a,b]和n 输出复合Simpson 积分值n S步1 .);()(;a x b f a f S nab h n ?-?-? 步2 对n k ,,2,1 =执行).(2;2);(4;2x f S S hx x x f S S h x x n n n n +?+?+?+?步3 n n S hS ??6步4 输出n SRomberg 积分法：根据已知数据对其进行多项式拟合得出p(x);f(x)?p(x); 输入被积函数f(x)，积分区间端点a,b,允许误差ε 输出 Romberg 积分值n R 2 步1 .0;0;0;0));()((2;1111?===+?-?k R C S b f a f hT a b h 步2 反复执行步3→步9. 步3 .2;0h a x S +步4 反复执行步5→步6. 步5 ;);(h x x x f S S +?+?步6 若x ≥b,则退出本层循环. 步7 执行.6316364;1511516;3134;2212212212212C C R S S C T T S S h T T -?-?-?+?步8 执行.1;;;;;2;2121212112+-?k k R R C C S S T T hh R R e 步9 若e ≤ε且k ≥5,则退出循环. 步10 .22R R n ? 步11 输出.2n R2、变步长梯形算法：功能求积分ba)(dx x f ，允许误差为ε。