西安交通大学_计算方法上机报告
计算方法上机报告_董晓壮
计算方法A上机报告学院（系）：电气工程学院学生姓名：陶然学号：授课老师：完成日期：2019年12月03日西安交通大学Xi'an Jiaotong University目录1 QR分解法求解线性方程组 (2)1.1 算法原理 (2)1.1.1 基于吉文斯变换的QR分解 (2)1.1.2 基于豪斯霍尔德变换的QR分解 (3)1.2 程序流程图 (4)1.2.1 基于吉文斯变换的QR分解流程图 (4)1.2.2 基于豪斯霍尔德变换的QR分解流程图 (5)1.3 程序使用说明 (5)1.3.1 基于吉文斯变换的QR分解程序说明 (5)1.3.2 基于豪斯霍尔德变换的QR分解程序说明 (7)1.4 算例计算结果 (8)2 共轭梯度法求解线性方程组 (10)2.1 算法原理 (10)2.2 程序流程图 (10)2.3 程序使用说明 (11)2.4 算例计算结果 (12)3 三次样条插值 (14)3.1 算法原理 (14)3.2 程序流程图 (16)3.3 程序使用说明 (17)3.4 算例计算结果 (19)4 龙贝格积分 (21)4.1 算法原理 (21)4.2 程序流程图 (22)4.3 程序使用说明 (23)4.4 算例计算结果 (24)结论 (26)1 QR 分解法求解线性方程组1.1 算法原理矩阵的QR 分解是指，可以将矩阵A 分解为一个正交矩阵Q 和一个上三角矩阵R 的乘积，实际中，QR 分解经常被用来解决线性最小二乘问题，分解情况如图1.1所示。
=⨯图1.1 QR 分解示意图本次上机学习主要进行了两个最基本的正交变换—吉文斯（Givens ）变换和豪斯霍尔德（Householder ）变换，并由此导出矩阵的QR 分解以及求解线性方程组的的方法。
1.1.1 基于吉文斯变换的QR 分解吉文斯矩阵也称初等旋转阵，如式（1.1）所示，它把n 阶单位矩阵I 的第,i j 行的对角元改为c ，将第i 行第j 列的元素改为s ，第j 行第i 列的元素改为s −后形成的矩阵。
(完整word版)计算方法A上机实验报告
计算方法A上机实验报告姓名：苏福班级：硕4020 学号：3114161019一、上机练习目的1）复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。
2）利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。
二、上机练习任务1）利用计算机语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。
2）掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。
3）写出上机练习报告。
三、上机题目1. 共轭梯度法求解线性方程组。
（第三章）2. 三次样条插值（第四章）3. 龙贝格积分（第六章）4. 四阶龙格－库塔法求解常微分方程的初值问题四、上机报告题目1：共轭梯度法求解线性方程组1．算法原理共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小值的问题。
从任意给定的初始点出发，沿一组关于矩阵A共轭的方向进行线性搜索，在无舍入误差的假定下，最多迭代n 次（其中n 为矩阵A 的阶数），就可求得二次函数的极小值，也就求得了线性方程组Ax b =的解。
定理：设A 是n 阶对称正定矩阵，则x *是方程组Ax b =的解得充分必要条件是x *是二次函数1()2TT f x x Ax b x =-的极小点，即 ()()min nx R Ax b f x f x **∈=⇔=共轭梯度法的计算公式：(0)(0)(0)()()()()(1)()()(1)(1)(1)()()()(1)(1)()k T k k k T k k k k k k k k T k k k T k k k k k d r b Ax r d d Ad xx d r b Ax r Ad d Ad d r d ααββ++++++⎧==-⎪⎪=⎪⎪=+⎪⎨=-⎪⎪⎪=-⎪⎪=+⎩2. 程序框图（1）编写共轭梯度法求解对称正定矩阵的线性方程组见附录（myge.m）：function x=myge(A,b)输入对称正定矩阵及对应的列向量，初始向量设为0，精度取为810 。
计算方法与实习上机实验报告
计算方法与实习上机实验报告一、引言本文旨在介绍和展示我们在“计算方法”课程中的实习上机实验环节所完成的一些关键任务和所取得的成果。
该实验课程的目标是让我们更深入地理解和应用各种计算方法，并在实际操作中提高我们的编程和问题解决能力。
二、实验内容与目标实验的主要内容是利用各种计算方法解决实际数学问题。
我们被要求使用编程语言（如Python或Java）来实现和解决这些问题。
这些问题包括使用牛顿法求解平方根，使用蒙特卡洛方法计算圆周率，以及使用最优化方法求解函数的最小值等。
实验的目标不仅是让我们掌握计算方法的基本理论，更是要让我们能够在实际操作中运用这些方法。
我们需要在实习过程中，通过与同伴们合作，共同解决问题，提高我们的团队合作能力和问题解决能力。
三、实验过程与问题解决策略在实验过程中，我们遇到了许多问题，如编程错误、理解困难和时间压力等。
我们通过相互讨论、查阅资料和寻求教师帮助等方式，成功地解决了这些问题。
例如，在实现牛顿法求解平方根时，我们一开始对导数的计算和理解出现了一些错误。
但我们通过查阅相关资料和讨论，最终理解了导数的正确计算方法，并成功地实现了牛顿法。
四、实验结果与结论通过这次实习上机实验，我们不仅深入理解了计算方法的基本理论，还在实际操作中提高了我们的编程和问题解决能力。
我们的成果包括编写出了能有效求解平方根、计算圆周率和求解函数最小值的程序。
这次实习上机实验非常成功。
我们的团队不仅在理论学习和实践操作上取得了显著的进步，还在团队合作和问题解决方面积累了宝贵的经验。
这次实验使我们对计算方法有了更深的理解和认识，也提高了我们的编程技能和解决问题的能力。
五、反思与展望回顾这次实验，我们意识到在实验过程中，我们需要更好地管理我们的时间和压力。
在解决问题时，我们需要更有效地利用我们的知识和资源。
在未来，我们希望能够更加熟练地运用计算方法，并能够更有效地解决问题。
我们也希望能够将所学的计算方法应用到更广泛的领域中，如数据分析、科学研究和工业生产等。
西安交通大学算法上机实验报告
《计算机算法设计与分析》上机实验报告姓名：班级：学号：日期：2016年12月23日算法实现题3-14 最少费用购物问题★问题描述：商店中每种商品都有标价。
例如，一朵花的价格是2元，一个花瓶的价格是5元。
为了吸引顾客，商店提供了一组优惠商品价。
优惠商品是把一种或多种商品分成一组，并降价销售。
例如，3朵花的价格不是6元而是5元。
2个花瓶加1朵花的优惠价格是10元。
试设计一个算法，计算出某一顾客所购商品应付的最少费用。
★算法设计：对于给定欲购商品的价格和数量，以及优惠价格，计算所购商品应付的最少费用。
★数据输入：由文件input.txt提供欲购商品数据。
文件的第1行中有1个整数B（0≤B≤5），表示所购商品种类数。
在接下来的B行中，每行有3个数C,K和P。
C表示商品的编码（每种商品有唯一编码），1≤C≤999；K表示购买该种商品总数，1≤K≤5；P是该种商品的正常单价（每件商品的价格），1≤P≤999。
请注意，一次最多可购买5*5=25件商品。
由文件offer.txt提供优惠商品价数据。
文件的第1行中有1个整数S（0≤S≤99），表示共有S种优惠商品组合。
接下来的S行，每行的第1个数描述优惠商品组合中商品的种类数j。
接着是j个数字对（C，K），其中C是商品编码，1≤C≤999；K表示该种商品在此组合中的数量，1≤K≤5。
每行最后一个数字P （1≤P≤9999）表示此商品组合的优惠价。
★结果输出：将计算出的所购商品应付的最少费用输出到文件output.txt。
输入文件示例输出文件示例Input.txt offer.txt output.txt2 2 147 3 2 1 7 3 58 2 5 2 7 1 8 2 10解：设cost（a，b，c，d，e）表示购买商品组合（a，b，c，d，e）需要的最少费用。
A[k]，B[k]，C[k]，D[k]，E[k]表示第k种优惠方案的商品组合。
offer （m）是第m种优惠方案的价格。
西安交大计算方法上机报告
从而得到计算的公式：
j 1, 2,..., n 1 j 1 j l i1 i 2,3,..., n i1 11 i 1 lik ukj j i , i 1,..., n, i 2,3,.., n ij ij k 1 i 1 1 l ( lkt ti ) k i 1,..., n, i 2,3,.., n ki ki t 1 ii
计算方法 上机实习题目报告
班级：材料 s3076 班 姓名：丁明帅 学号：3113305029
1：计算 S
100000
k 1
1 ,要求误差小于 106 ,给出实现算法。 k2
【实现思路】
设当 k 值为 i 时 S-Si<10-6,则只需要
S Si
100000 1
k i
1
2
1 l32 ln 2 1 ln 3 1
4 / 46
12 22
1n 2 n
nn
因此有
1 j 2 j 0 jj 0 0
ij li1
li 2
li ,i -1
根据定义知道 L 矩阵的斜对角线上的值都为 1，且 L 矩阵的第一行与原矩阵 A 的第一行 相同，因此可以根据公式先计算 U 的第一行，然后计算 L 的第一列；以后的第 i 步先计算 U 的第 i 行， 然后计算 L 的第 i 列 （U 的第 n 行不作计算） 。 然后把最后的 U 的第 n 行计算出来。
【算法依据】
列主元高斯消元法的过程可以将方程组系数简化为系数矩阵与 b 矩阵， 从而利用方程组 对系数扩展矩阵进行消元。 在消元的过程中矩阵的行向量之间可以变换， 但列向量不能变化。 在进行压缩矩阵的求解中还需要人为的将因调整行向量所导致的列向量的变化调整回来。 Matlab 对于矩阵的处理非常容易，结合 for 循环语句可以实验本题大规模方程组的求解。
计算方法与实习上机报告
计算方法与实习——上机报告学院：电子工程学院学号：姓名：刘波2015.1.4计算方法与实习上机报告习题一：1 舍入误差及稳定性一、实验目的（1）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；（2）通过上机计算，了解舍入误差所引起的数值不稳定性二、实验内容1、用两种不同的顺序计算1000021n n -=∑，分析其误差的变化 2、已知连分数()101223//(.../)n n a f b b a b a a b =++++，利用下面的算法计算f ： 11,i n n i i i a d b d b d ++==+(1,2,...,0)i n n =-- 0f d = 写一程序，读入011,,,...,,,...,,n n n b b b a a 计算并打印f 3、给出一个有效的算法和一个无效的算法计算积分1041nn x y dx x =+⎰ (0,1,...,10)n = 4、设2211N N j S j ==-∑，已知其精确值为1311221N N ⎛⎫-- ⎪+⎝⎭ （1）编制按从大到小的顺序计算N S 的程序 （2）编制按从小到大的顺序计算N S 的程序（3）按两种顺序分别计算10001000030000,,,S S S 并指出有效位数三、实验步骤、程序设计、实验结果及分析1、用两种不同的顺序计算1000021n n -=∑，分析其误差的变化 （1）实验步骤：分别从1~10000和从10000~1两种顺序进行计算，应包含的头文件有stdio.h 和math.h（2）程序设计：a.顺序计算#include<stdio.h>#include<math.h>void main(){double sum=0;int n=1;while(1){sum=sum+(1/pow(n,2));if(n%1000==0)printf("sun[%d]=%-30f",n,sum);if(n>=10000)break;n++;}printf("sum[%d]=%f\n",n,sum);}b.逆序计算#include<stdio.h>#include<math.h>void main(){double sum=0;int n=10000;while(n!=0){sum=sum+(1/pow(n,2));if(n%200==0)printf("sum[%d]=%-10f",n,sum);if(n<1)break;n--;}printf("sum[%d]=%f\n",n,sum);}（3）实验结果及分析：程序运行结果：a.顺序计算b.逆序计算结果分析：两种不同顺序计算结果是一样的，顺序计算误差从一开始就很小，而逆序计算误差最开始十分大，后来结果正确。
西安交大 计算方法B上机作业
计算方法（B ）上机作业一、三次样条拟合某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:米)如下表所示：(1)(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图； 解：1、算法实现的思想及依据题目（1）为曲线拟合问题多项式插值、分段插值和最小二乘法。
多项式插值，随着插值数据点的数目增多，误差也会随之增大，因此不选用。
最小二乘法适于数据点较多的场合，在此也不适用。
故选用分段插值。
分段插值又分为分段线性插值、分段二次插值、三次样条插值及更高阶的多项式插值。
由本题的物理背景知，光缆正常工作时各点应该是平滑过渡，因此至少选用三次样条插值法。
对于更高阶的多项式插值，由于“龙格现象”而不选用。
题目（2）求光缆长度，即求拟合曲线在0到20的长度，对弧长进行积分即可。
光缆长度的第一型线积分表达式为190k kk l +==∑⎰。
2、算法实现的结构参考教材给出的SPLINEM 算法和TTS 算法，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出来并赋值给右端向量d ，再根据TSS 解法求解M 。
光缆长度的第一型线积分表达式为190k kk l +==∑⎰。
3、程序运行结果及分析图1.1三种边界条件下三次样条插值图1.2光缆长度4、MATLAB代码：1）自己编程实现代码clear;clc;I=input('你想使用第几种边界条件？请输入1、2、3之一: ');x=0:20;y=[9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.22 9.15 7.90 7.95 8.86 9.81 10.8 10.93];plot(x,-y,'k.','markersize',15)%y为深度，取负号hold on%% 计算一阶差商y1=ones(1,21);for i=2:1:21y1(i)=(y(i)-y(i-1))/(x(i)-x(i-1));end%% 计算二阶差商y2=ones(1,21);for i=3:1:21y2(i)=(y1(i)-y1(i-1))/(x(i)-x(i-2));end%% 计算三阶差商y3=ones(1,21);for i=4:1:21y3(i)=(y2(i)-y2(i-1))/(x(i)-x(i-3));end%% 选择边界条件（I）if I==1d(1)=0;d(21)=0;a(21)=0;c(1)=0;% 第一个点和最后一个点的二阶差商为0 endif I==2d(1)=6*y1(1);d(21)=-6*y1(21);a(1)=1;c(1)=1;endif I==3d(1)=-12*y3(1);d(21)=-12*y3(21);a(21)=-2;c(1)=-2;%endfor i=2:20d(i)=6*y2(i+1);end%% 构造带状矩阵求解（追赶法）b=2*ones(1,21);a=0.5*ones(1,21);%a(21)=-2;c=0.5*ones(1,21);%c(1)=-2;u(1)=b(1);r(1)=c(1);%% 追yz(1)=d(1);for i=2:21l(i)=a(i)/u(i-1);u(i)=b(i)-l(i)*r(i-1);r(i)=c(i);yz(i)=d(i)-l(i)*yz(i-1);end%% 赶xg(21)=yz(21)/u(21);for i=20:-1:1xg(i)=(yz(i)-r(i)*xg(i+1))/u(i);endM=xg;%%所有点的二阶导数值%% 求函数表达式并积分t=1;h=1;N=1000x1=0:20/(N-1):20;length=0;for i=1:Nfor j=2:20if x1(i)<=x(j)t=j;break;elset=j+1;endendf1=x(t)-x1(i);f2=x1(i)-x(t-1);S(i)=(M(t-1)*f1^3/6/h+M(t)*f2^3/6/h+(y(t-1)-M(t-1)*h^2/6)*f1+(y(t)-M(t)*h^2/6)* f2)/h;Sp(i)=-M(t-1)*f1^2/2/h+M(t)*f2^2/2/h+(y(t)-y(t-1))/h-(M(t)-M(t-1))*h/6;length(i+1)=sqrt(1+Sp(i)^2)*(20/(N-1))+length(i);%第一类线积分endfigure(1);plot(x1,-S,'r-')%深度曲线griddisp(['第',num2str(I),'种边界条件下长度',num2str(length(N+1)),'米'])axis fill;xlabel('测点/米');ylabel('深度/米');title('三次样条曲线拟合');legend('数据点','拟合曲线',3);二、最小二乘近似假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。
西安交通大学第二次汇编上机实验报告
西安交通大学第二次汇编上机实验报告
实验时间5.31姓名
学号
实验内容概述主要对第五章所学的循环了分支结构进行了实践。
第一题循环结构比大小，把最大值的最大值的相对位置存储在数组之后的两位第二题循环结构找两数组之间的相同数字第三题分支结构对数组里出现的元素进行计数实验步骤
一（1）对数据求源码（负数取反后加一）（2）取每个数的绝对值（3）从前向后两两比较，如果当前数字比max的记录值大，那么更新当前数字的max，并存储当前位置pos二（1）根据数组a的长度设置外循环次数（2）根据数组b 的长度设置外循环次数（3）外循环遍历a的每一个元素，将a中每一个元素与b中元素比较，如果相同，则终止当前内循环，将该值存储在c数组中三（1）定义一个分支结构，先判断当前值，如果相等，则跳转至对应子程序：使存储该数字出现次数的位置上加一（2）程序运行结束时，在dos窗口输出每个元素出现的次数，
实验中遇到的问题及其解决方法1、mov指令用错当时报错后，查找了mov 指令的课本内容，想起来存储单元之间不能直接赋值，已解决2、逻辑正确，数据错误当时确认逻辑无误之后，我的内心非常纠结，甚至出现了想要重装软件的想法，后来利用u和t指令，一步步找错，发现是比例因子出现错误，将inccx 改为addcx，2后解决实验结果
结果均正确
2021
年
5月
31日
建议无备注无。
西安交通大学计算方法A上机作业
计算方法（A）大作业姓名：班级：专业：学号：共轭梯度法一、算法原理共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小化的问题，因此从任意给定的初始点出发，沿一组关于矩阵A的共轭方向进行线性搜索，在无舍入无差的假定下，最多迭代n(其中n为矩阵A的阶数）次就可求得二次函数的极小点，也就求得了线性方程组Ax=B的解。
下述定理给出了求系数矩阵A是对称正定矩阵的线性方程组Ax=b的解与求二次函数f(x)=12x T Ax−b T x极小点的等价性。
定理3.4.1设A是n阶对称正定矩阵，则x∗是方程组Ax=b的解的充分必要条件是x∗是二次函数f(x)=12x T Ax−b T x的极小点，即Ax∗=b⟺f(x∗)=minx∈R nf(x)证明：充分性.设x∗是f(x)的极小点，则由无约束最优化问题最优解的必要条件知∇f(x∗)=Ax∗−b即x∗是方程组Ax=b的解。
必要性. 若x∗是方程组Ax=b的解，即A x∗=b，注意到A是对称正定矩阵，故∀x∈R n有f(x)−f(x∗)=12x T Ax−b T x−12x T Ax∗+b T x∗=12(x T Ax−2b T x+x∗T Ax∗)−x∗T Ax∗+b T x∗=12(x T Ax−2(Ax∗)T x+x∗T Ax∗)−(Ax∗−b)T x∗=12(x−x∗)T A(x−x∗)≥0即x∗是f(x)的极小点，进而由A是正定矩阵知，x∗是f(x)的最小点。
证毕。
共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式：x(k+1)=x(k)+αk d(k)产生的迭代序列x(1),x(2),x(3)…在无舍入误差假定下，最多经过n次迭代，就可求得f(x) 的最小值，也就是方程Ax=b的解。
共轭梯度法中关键的两点是，确定迭代格式中的搜索向量d(k)和最佳步长αk （αk≥0）。
实际上，搜索方向d(k)是关于矩阵A的共轭向量，在迭代中逐步构造之。
西安交大计算方法B2017大作业
计算方法B上机报告姓名：学号：班级：学院：任课教师：2017年12月29日题目一：1.1题目内容某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:(1)(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；1.2 实现题目的思想及算法依据首先在题目（1）中要实现的是数据的拟合，显然用到的是我们在第三章中数据近似的知识内容。
多项式插值时，这里有21个数据点，则是一个20次的多项式，但是多项式插值随着数据点的增多，会导致误差也会随之增大，插值结果会出现龙格现象，所以不适用于该题目中点数较多的情况。
为了避免结果出现大的误差，同时又希望尽可能多地使用所提供的数据点，提高数据点的有效使用率，这里选择分段插值方法进行数据拟合。
分段插值又可分为分段线性插值、分段二次插值和三次样条插值。
由于题目中所求光缆的现实意义，而前两者在节点处的光滑性较差，因此在这里选择使用三次样条插值。
根据课本SPLINEM 算法和TSS 算法，采用第三种真正的自然边界条件，在选定边界条件和选定插值点等距分布后，可以先将数据点的二阶差商求出并赋值给右端向量d ，再根据TSS 解法求解三对角线线性方程组从而解得M 值。
求出M 后，对区间进行加密，计算200个点以便于绘图以及光缆长度计算。
对于问题（2），使用以下的公式:20=()L f x ds ⎰20(f x =⎰191(k kk f x +==∑⎰1.3 算法结构1. For n i ,,2,1,0⋅⋅⋅=1.1 i i M y ⇒2. For 2,1=k2.1 For k n n i ,,1, -=2.1.1 i k i i i i M x x M M ⇒----)/()(13. 101h x x ⇒-4. For 1-,,2,1n i =4.1 11++⇒-i i i h x x4.2 b a c c h h h i i i i i i ⇒⇒-⇒+++2;1;)/(11 4.3 i i d M ⇒+165. 0000;;c M d M d n n ⇒⇒⇒λn n n b a b ⇒⇒⇒2;;20μ6. 1111,γμ⇒⇒d b7. For m k ,,3,2 = ! 获取M 的矩阵元素个数，存入m7.1 k k k l a ⇒-1/μ 7.2 k k k k c l b μ⇒⋅-1- 7.3 k k k k l d γγ⇒⋅-1- 8. m m m M ⇒μγ/9. For 1,,2,1 --=m m k9.1 k k k k k M M c ⇒⋅-+μγ/)(110. k ⇒1 ! 获取x 的元素个数存入s 11. For 1,,2,1-=s i11.1 if i x x ≤~then k i ⇒；break else k i ⇒+112. xx x x x x h x x k k k k ˆ~;~;11⇒-⇒-⇒--- y h x h M y x h M y x M x M k k k k k k ~/]ˆ)6()6(6ˆ6[2211331⇒-+-++---1.4 matlab 源程序n=20; x=0:n;y=[9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.22 9.15 7.90 7.95 8.86 9.81 10.80 10.93];M=y; %用于存放差商，此时为零阶差商 h=zeros(1,n+1); c=zeros(1,n+1); d=zeros(1,n+1); a=zeros(1,n+1); b=2*ones(1,n+1); h(2)=x(2)-x(1);for i=2:n %书本110页算法SPLINEM h(i+1)=x(i+1)-x(i);c(i)=h(i+1)/(h(i)+h(i+1)); a(i)=1-c(i); enda(n+1)=-2; %计算边界条件c(0),a(n+1),采用的是第三类边界条件 c(1)=-2;for k=1:3 %计算k 阶差商 for i=n+1:-1:k+1M(i)=(M(i)-M(i-1))/(x(i)-x(i-k)); endif(k==2) %计算2阶差商 d(2:n)=6*M(3:n+1); %给d 赋值 endif(k==3)d(1)=(-12)*h(2)*M(4); %计算边界条件d(0),d(n),采用的是第三类边界条件 d(n+1)=12*h(n+1)*M(n+1); end endl=zeros(1,n+1); r=zeros(1,n+1); u=zeros(1,n+1); q=zeros(1,n+1); u(1)=b(1); r(1)=c(1); q(1)=d(1);for k=2:n+1 %利用书本49页算法TSS求解三对角线性方程组r(k)=c(k);l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*r(k-1);q(k)=d(k)-l(k)*q(k-1);endp(n+1)=q(n+1)/u(n+1);for k=n:-1:1p(k)=(q(k)-r(k)*p(k+1))/u(k);endfprintf('三对角线性方程组的解为：');disp(p);%求拟合曲线x1=0:0.1:20; %首先对区间进行加密,增加插值点n1=10*n;x2=zeros(1,n1+1);x3=zeros(1,n1+1);s=zeros(1,n1+1);for i=1:n1+1for j=1:nif x1(i)>=x(j)&&x1(i)<=x(j+1) %利用书本111页算法EVASPLINE求解拟合曲线s(x)h(j+1)=x(j+1)-x(j);x2(i)=x(j+1)-x1(i);x3(i)=x1(i)-x(j);s(i)=(p(j).*(x2(i)).^3/6+p(j+1).*(x3(i)).^3/6+(y(j)-p(j).*((h(j+1)).^2/6)).*x2( i)+...(y(j+1)-p(j+1).*(h(j+1)).^2/6).*x3(i))/h(j+1);endendendplot(x,-y,'x') %画出插值点hold onplot(x1,-s) %画出三次样条插值拟合曲线hold ontitle('三次样条插值法拟合电缆曲线');xlabel('河流宽度/m');ylabel('河流深度/m');Length=0;for i=1:n1L=sqrt((x1(i+1)-x1(i))^2+(s(i+1)-s(i))^2); %计算电缆长度Length=Length+L;endfprintf('电缆长度(m)=');disp(Length);1.5 结果与说明铺设海底光缆的曲线如图1.1所示图1. 1三次样条插值法拟合海底光缆曲线由上图可以看出，所得到的曲线光滑，能够较好得反映实际的河沟底部地势形貌。
计算方法上机实验报告
《计算方法》上机实验报告班级:XXXXXX小组成员：XXXXXXXXXXXXXXXXXXXXXXXXXXXX任课教师：XXX二〇一八年五月二十五日前言通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB 程序的编写.以下为本次上机实验报告，按照实验内容共分为六部分.实验一：一、实验名称及题目： Newton 迭代法例2。
7（P38):应用Newton 迭代法求在附近的数值解，并使其满足.二、解题思路：设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交点的横坐标)(')(0001x f x f x x -=,称1x 为'x 的一次近似值，过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标)(')(1112x f x f x x -=称2x 为'x 的二次近似值，重复以上过程，得'x 的近似值序列{}n x ,把)(')(1n n n n x f x f x x -=+称为'x 的1+n 次近似值，这种求解方法就是牛顿迭代法。
三、Matlab 程序代码：function newton_iteration （x0,tol) syms z ％定义自变量 format long ％定义精度 f=z*z*z —z-1;f1=diff(f ）；%求导 y=subs(f ，z ，x0);y1=subs （f1,z,x0)；%向函数中代值 x1=x0-y/y1; k=1;while abs(x1-x0）>=tol x0=x1；y=subs(f,z,x0)； y1=subs （f1,z,x0）; x1=x0—y/y1;k=k+1; endx=double(x1) K四、运行结果：实验二：一、实验名称及题目：Jacobi 迭代法例3。
西交计算方法A上机大作业
计算方法A 上机大作业1. 共轭梯度法求解线性方程组算法原理：由定理3.4.1可知系数矩阵A 是对称正定矩阵的线性方程组Ax=b 的解与求解二次函数1()2TT f x x Ax b x =- 极小点具有等价性，所以可以利用共轭梯度法求解1()2TT f x x Ax b x =-的极小点来达到求解Ax=b 的目的。
共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式：(1)()()k k k k x x d α+=+产生的迭代序列(1)(2)(3)x x x ，，，... 在无舍入误差假定下，最多经过n 次迭代，就可求得()f x 的最小值，也就是方程Ax=b 的解。
首先导出最佳步长k α的计算式。
假设迭代点()k x 和搜索方向()k d 已经给定，便可以通过()()()()k k f x d φαα=+的极小化()()min ()()k k f x d φαα=+来求得，根据多元复合函数的求导法则得:()()()'()()k k T k f x d d φαα=∇+令'()0φα=，得到:()()()()k T k k k T k r d d Adα= ，其中()()k k r b Ax =-然后确定搜索方向()k d 。
给定初始向量(0)x 后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向(0)(0)(0)(0)()dr f x b Ax ==-∇=- 。
令(1)(0)00x x d α=+其中(0)(0)0(0)(0)T T r d d Adα=。
第二次迭代时，从(1)x 出发的搜索方向不再取(1)r ，而是选取(1)(1)(0)0dr d β=+，使得(1)d 与(0)d 是关于矩阵A 的共轭向量，由此可求得参数0β:(1)(0)0(0)(0)T T r Ad d Adβ=-然后从(1)x 出发，沿(1)d 进行搜索得到(2)(1)(1)1x x d α=+设已经求出(1)()()k k k k x x d α+=+，计算(1)(1)k k r b Ax ++=-。
西安交通大学 计算方法上机报告
计算方法上机报告姓名：学号：班级：机械硕4002上课班级：02班说明：本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。
1. 对以下和式计算：∑∞⎪⎭⎫ ⎝⎛+-+-+-+=0681581482184161n n n n S n，要求：① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算；（1） 算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 142111416818485861681n n n a n n n n n ε⎛⎫=---<< ⎪+++++⎝⎭； 2、为了保证计算结果的准确性，写程序时，从后向前计算； 3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)（2）算法结构1. ;0=s⎪⎭⎫⎝⎛+-+-+-+=681581482184161n n n n t n； 2. for 0,1,2,,n i =⋅⋅⋅ if 10m t -≤end;3. for ,1,2,,0n i i i =--⋅⋅⋅;s s t =+（3）Matlab源程序clear; %清除工作空间变量clc; %清除命令窗口命令m=input('请输入有效数字的位数m='); %输入有效数字的位数s=0;for n=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));if t<=10^(-m) %判断通项与精度的关系break;endend;fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值for i=n-1:-1:0t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));s=s+t; %求和运算ends=vpa(s,m) %控制s的精度（4）结果与分析当保留11位有效数字时，需要将n值加到n=7，s =3.1415926536；当保留30位有效数字时，需要将n值加到n=22,s =3.14159265358979323846264338328。
西安交通大学计算方法上机作业
计算方法上机作业1.对以下和式计算：0142118184858616n n S n n n n ∞=⎛⎫=--- ⎪++++⎝⎭∑，要求： (1)若只需保留11个有效数字，该如何进行计算； (2)若要保留30个有效数字，则又将如何进行计算；（1）解题思想和算法实现：根据保留有效位数的要求，可以由公式得出计算精度要求。
只需要很少内存，时间复杂度和d 呈线性，不需要高浮点支持。
先根据while 语句求出符合精度要求的n 值的大小，然后利用for 语句对这n 项进行求和，输出计算结果及n 值大小即可。
（2）matlab 源程序：保留11位有效数字时； clear clcformat long n=0;sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); while sum>=5*10^(-11); n=n+1;sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); endfor i=0:n-1;sum=sum+1/(16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); endvpa(sum,11) n保留30位有效数字时； clear clcformat long n=0;sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); while sum>=5*10^(-30); n=n+1;sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); endfor i=0:n-1;sum=sum+1/(16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); endvpa(sum,30) n（3）实验结果分析图1.1 保留11位有效数字的n值及计算结果图图1.2 保留30位有效数字的n值及计算结果图由计算结果可知，通过合理的误差控制，分别通过7次和22次循环，可以实现题目所要求的精确度。
西安交大计算方法
西安交通大学计算方法上机实验班级：(xxx)姓名：(xxx)学号：21116010041.按两种顺序计算y，哪个接近真值？Y = 1000 + + + … +用java 语言编写：public class Add {public static void main(String[] args){double s=0,y=1000;for(double a=1001.0;a<=2000.0;a++){y+=1.0/a;}for(double a=2000.0;a>=1001.0;a--){s+=1.0/a;}s=s+1000;System.out.println("正序和"+s);System.out.println("逆序和"+y);}}运行结果：结论：显然假设是double类型的数据时，先算大数的过程吃掉了末尾的小数被进位所埋没，导致了大数吃小数的误差，按从小到大（从右向左）的计算顺序所得的结果与真值相近，而按从大到小（从左到右）的计算顺序所得的结果与真值的误差较大。
1-18.设(x) = 1 + x + + + … + , 计算(-5)和1/(5)，哪个接近？解法一：用JAVA 语言编写：public class second{ public static void main(String[] args){double s1=1 ,s2=1;double e=1,sum=1; //e的初值为1，sum用来存放n！int a=1;while(sum<Math.pow(10, 1000000)){sum=a*sum;e=1.0/sum+e;a++;}double b=1.0/(e*e*e*e*e);System.out.println("较为精确的值1/e^5="+b);for(int i=1;i<=24;i++){s1+=cimi1(i);s2+=cimi2(i);}s1=1.0/s1;System.out.println("1/S24(5)="+s1);System.out.println("S24(-5)="+s2);}public static double cimi1(int ai){double xi=1;for(int i=ai;i>=1;i--){xi=xi*(5.0/i);}return xi;}public static double cimi2(int ai){double xi=1;for(int i=ai;i>=1;i--){xi=xi*(-5.0/i);}return xi;}}运行结果：解法二：用matlab编程并运行，如下：(1)计算(-5)运行结果如下：(2)计算1/(5)运行结果如下：而的真是结果为0.006737946比较得1/(5)的计算结果与真实值更接近解法三：也可以用C++编写：#include "stdafx.h"#include"stdio.h"#include "iostream"using namespace std;int main(int argc, char* argv[]){ int func1(int );double func2(int);double y=0;int i;for(i=1;i<25;i++){ int z=func1(i);double e=func2(i);y+=z/e;}cout<<"----------------------------------------"<<endl;cout<<"1/S(5)的运算结果是："<<" "<<1.0/(y+1)<<endl;cout<<"----------------------------------------"<<endl;return 0;}int func1(int x){int y=1;int k;for (k=0;k<x;k++)y*=5;return y;}double func2(int n){double y=1;int j;for (j=1;j<=n;j++)y*=j;return y;}运行结果如下图：结论：通过比较上述的几种编程结果，可以看出1/S(5),更接近真实值，而且用matlab更为简便，可以直接利用函数库，并可以轻松的嵌入秦九韶算法，大大减少运算量和时间。
西安交大计算方法b大作业
西安交大计算方法b大作业《计算方法B》上机实验报告学院：班级: 姓名：学号：机械工程学院2021年12月22日11.计算以下和式：S??1nn?016?211??4????8n?18n?48n?58n?6?，要求： ??(1)若保留11个有效数字，给出计算结果，并评价计算的算法；(2)若要保留30个有效数字，则又将如何进行计算。
实现思想：以上问题出现了近似数相减的问题，为了减小误差，可分别求得减数之和以及被减数之和，最后将两者相减。
另外，减数与被减数求和均为同号计算，按照绝对值递增顺序相加可减小舍入误差。
此题中对有效数字有要求，因而计算时首先需要根据有效数字位数计算得出迭代次数，以保证计算值的精度。
源程序：m=input('输入有效数字个数m=');s0=1;s1=0;s2=0;n=0; %判断迭代次数while s0>=0.5*10^-(m-1)s0=4/(16^n*(8*n+1))-2/(16^n*(8*n+4))-1/(16^n*(8*n+5))-1/(16^n*(8*n+6)); n=n+1; end%分别求解各项并求和 for k=n-1:-1:0a1=4/(16^k*(8*k+1)); a2=2/(16^k*(8*k+4)); a3=1/(16^k*(8*k+5)); a4=1/(16^k*(8*k+6)); s1=a1+s1; s2=a4+a3+a2+s2; endS=vpa(s1-s2,m)2实验结果：11位有效数字计算结果如图1所示；30为有效数字计算结果如图2所示。
图1.11位有效数字计算结果图2.30为有效数字计算结果31. 某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:米)如下表所示：分点深度分点深度分点深度 0 9.01 7 11.18 14 9.15 1 8.96 8 12.26 157.90 2 7.96 9 13.28 16 7.95 3 7.97 10 13.32 17 8.86 4 8.02 11 12.61 18 9.81 5 9.05 12 11.29 19 10.80 6 10.13 13 10.22 20 10.93 (1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；算法思想：由于题中所给点数为20，若采用高次多项式插值将产生很大的误差，所以拉格朗日或牛顿并不适用。
计算方法上机报告
计算方法上机报告硕5028李兰鑫31150350051. 计算以下和式：0142118184858616nn S n n n n ∞=⎛⎫=--- ⎪++++⎝⎭∑，要求： (1)若保留11个有效数字，给出计算结果，并评价计算的算法； (2)若要保留30个有效数字，则又将如何进行计算。
1.1算法思想由题可得S<(1/16^n)*(4/(8*n+1))保留11个有效数字，即误差ε<0.5*10^（1-11） 保留30个有效数字，即误差ε<0.5*10^（1-30） 1.2源程序%有效数字为11的情况digita=11;（第二问中改为30即可） for i=1:1:1000S0=(1/(16^(i-1)))*(4/(8*i-7)); S1=(1/16^i)*(4/(8*i+1)); if ((S0-S1)<0.5*10^(1-digita)) n=i; break ; end end S=0;digits(digita);for i=0:1:nS=S+vpa(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); end1.3运行结果2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测(1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；2.1算法思想由题可知应使用多项式插值，考虑到本例中的拟合数据点较多，使用多项式插值时会因龙格现象造成误差较大，故决定采用三次样条插值。
2.2源程序%------三次样条插值多项式------clc;clear;%-----插值函数和插值点-------n=20;%插值点个数11i=0:n;%子元素从0到nx=10*i;f=[9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.229.15 7.90 7.95 8.86 9.81 10.80 10.93];%--输出插值点---fprintf('当n=%d时\n',n)fprintf('插值点x坐标为：\n')disp(x);fprintf('插值点y坐标f（x）为：\n')disp(f);%----计算差商表，并将有用的差商进行选择和保存-----y=f;%存放差商d=zeros(1,n+1);%存放三弯矩方程组的右端项for i=1:n%计算第i阶差商for j=n+1:-1:i+1y(j)=(y(j)-y(j-1))/(x(j)-x(j-i));endif i==2%将第2阶差商赋值给dd(2:n)=6.*y(3:n+1);endif i==3%三次样条插值的第三类边界条件d(1)=-12*(x(2)-x(1))*y(i+1);d(n+1)=12*(x(n+1)-x(n))*y(n+1); endend%-------用追赶法求三弯矩方程组------a=ones(1,n+1);%存储下三角的非零元素b=2*ones(1,n+1);%存储对角元素c=ones(1,n+1);%存储上三角的非零元素for i=2:na(i)=(x(i)-x(i-1))/(x(i+1)-x(i-1));c(i)=1-a(i);enda(n+1)=-2;c(1)=-2;%-------计算矩阵L和U-----------%初始化L和Ul=zeros(1,n+1);u=zeros(1,n+1);r=zeros(1,n+1);u(1)=b(1);r(1)=c(1);for k=2:n+1l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*r(k-1);r(k)=c(k);end%--------使用追赶法进行求解----------M=zeros(1,n+1);%方程组的解M(1)=d(1);for k=2:n+1M(k)=d(k)-l(k)*M(k-1);endM(n+1)=M(n+1)/u(n+1);for k=n:-1:1M(k)=(M(k)-r(k)*M(k+1))/u(k);endfprintf('三弯矩方程组的解M为：\n')disp(M);%------计算f（x）、N(x)和S(x)-------Sy1=zeros(1,200);x1=zeros(1,200);for i=1:200x1(i)=i-1;endfor i=1:200for j=1:nif x1(i)>=x(j)&&x1(i)<=x(j+1)temp=x(j+1)-x(j);Sy1(i)=M(j).*(x(j+1)-x1(i)).^3/6/temp+M(j+1).*(-x(j)+x1(i)).^3/6/temp+(f(j)-M(j)*te mp^2/6).*(x(j+1)-x1(i))/temp+(f(j+1)-M(j+1)*temp^2/6).*(-x(j)+x1(i))/temp;endendendfprintf('由三次样条插值方法计算的各年的产量为：\n')disp(Sy1)plot(x,f,'*')%画出由插值点构成的曲线hold onplot(x1,Sy1,'r')%用红线画出三次样条插值拟合的曲线hold onlegend('三次样条插值');%------求解电缆长度------Cordinate=[0:20/(length(Sy1)-1):20;Sy1];diff_Cordinate=[diff(Cordinate(1,:));diff(Cordinate(2,:))];disp('电缆长度为：');anwser=sum(sqrt(diff_Cordinate(1,:).*diff_Cordinate(1,:)+diff(Cordinate(2,:)).*diff( Cordinate(2,:))))2.3运行结果由图可知，三次样条插值拟合数据点比较好，没有龙格现象。
计算方法实验上机报告(完整版)
简单迭代法#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) pow(12*x+sin(x)-1,1.0/3)void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果牛顿迭代法一#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) x-(pow(x,3)-sin(x)-12*x+1)/(3*pow(x,2)-cos(x)-12) void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果埃特金加速法#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) (pow(x,3)-sin(x)+1)/12void main(){int i;double x1=x0,x2=x0;double z,y;printf("k\tx1\tx2\txk\n");for(i=0;i<MAXREPT;i++){if(i==0)printf("%d\t\t\t%g\n",i,x2);elseprintf("%d\t%g\t%g\t%g\n",i,y,z,x2);y=G(x1);z=G(y);x2=z-((z-y)*(z-y))/(z-2*y+x1);if (fabs(x2-x1)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x2,i);return;}x1=x2;}printf("AFTER %d repeate,no solved.\n",MAXREPT);} 结果牛顿迭代法二#include<stdio.h>#include<math.h>#define x0 1.5#define MAXREPT 1000#define EPS 1E-6#define G(x) x-(pow(x,3)+pow(x,2)-3*x-3)/(3*pow(x,2)+2*x-3) void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果弦截法#include<stdio.h>#include<math.h>#define x0 0#define x1 1.5#define MAXREPT 1000#define EPS 1E-6#define G(x) pow(x,3)+pow(x,2)-3*x-3void main(){int i;double x_k=x0,x_k1=x1,x_k2=0;double y,z;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k2);y=G(x_k);z=G(x_k1);x_k2=x_k1-(z*(x_k1-x_k))/(z-y);if (fabs(x_k2-x_k1)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k2,i);return;}x_k=x_k1;x_k1=x_k2;}printf("AFTER %d repeate,no solved.\n",MAXREPT); } 结果高斯顺序消元法#include<stdio.h>#include<math.h>#define N 4static double aa[N][N+1]={{2,4,0,1,1},{3,8,2,2,3},{1,3,3,0,6},{2,5,2,2,3}}; int gauss(double a[][N+2],double x[]);void putout(double a[][N+2]);void main(){int i,j,det;double a[N+1][N+2],x[N+1];for(i=1;i<=N;i++)for(j=1;j<=N+1;j++)a[i][j]=aa[i-1][j-1];det=gauss(a,x);if(det!=0)for(i=1;i<=N;i++)printf(" x[%d]=%g",i,x[i]);printf("\n");}int gauss(double a[][N+2],double x[]){int i,j,k;double c;putout(a);for(k=1;k<=N-1;k++){ if(fabs(a[k][k])<1e-17){printf("\n pivot element is 0.fail!\n");return 0;}for(i=k+1;i<=N;i++){c=a[i][k]/a[k][k];for(j=k;j<=N+1;j++){a[i][j]=a[i][j]-c*a[k][j];}}putout(a);}if(fabs(a[N][N])<1e-17){printf("\n pivot element is 0.fail!\n");return 0;}for(k=N;k>=1;k--){x[k]=a[k][N+1];for(j=k+1;j<=N;j++){x[k]=x[k]-a[k][j]*x[j];}x[k]=x[k]/a[k][k];}return(1);}void putout(double a[][N+2]){for(int i=1;i<=N;i++){for(int j=1;j<=N+1;j++)printf("%-15g",a[i][j]);printf("\n");}printf("\n");}结果雅克比迭代法#include<stdio.h>#include<math.h>#define N 5#define EPS 0.5e-4static double aa[N][N]={{4,-1,0,-1,0},{-1,4,-1,0,-1},{0,-1,4,-1,0},{-1,0,-1,4,-1},{0,-1,0,-1,4}}; static double bb[N]={2,1,2,1,2};void main(){int i,j,k,NO;double a[N+1][N+1],b[N+1],x[N+1],y[N+1];double d,sum,max;for(i=1;i<=N;i++){for(j=1;j<=N;j++)a[i][j]=aa[i-1][j-1];b[i]=bb[i-1];}printf("\n 请输入最大迭代次数（尽量取大值！）NO:");scanf("%d",&NO);printf("\n");for(i=1;i<=N;i++)x[i]=0;k=0;printf(" k",' ');for(i=1;i<=N;i++)printf("%8cx[%d]",' ',i);printf("\n 0");for(i=1;i<=N;i++)printf("%12.8g",x[i]);printf("\n");do{for(i=1;i<=N;i++){sum=0.0;for(j=1;j<=N;j++)if(j!=i) sum=sum+a[i][j]*x[j];y[i]=(-sum+b[i])/a[i][i];}max=0.0;for(i=0;i<=N;i++){d=fabs(y[i]-x[i]);if(max<d) max=d;x[i]=y[i];}printf("%6d",k+1);for(i=1;i<=N;i++)printf("%12.8g",x[i]);printf("\n");k++;}while((max>=EPS)&&(k<NO));printf("\nk=%d\n",k);if(k>=NO) printf("\nfail!\n");elsefor(i=1;i<=N;i++)printf("x[%d]=%g\t",i,x[i]);}结果拉格朗日插值多项式#include<stdio.h>#include<math.h>#define N 4doublex[N]={0.56160,0.56280,0.56401,0.56521},y[N]={0.82741,0.82659,0.82577,0.82495}; void main(){double x=0.5635;double L(double xx);double lagBasis(int k,double xx);void output();output();printf("\n近似值L(%g)=%g\n",x,L(x));}double lagBasis(int k,double xx){double lb=1;int i;for(i=0;i<N;i++)if(i!=k) lb*=(xx-x[i])/(x[k]-x[i]);return lb;}double L(double xx){double s=0;int i;for(i=0;i<=N;i++)s+=lagBasis(i,xx)*y[i];return s;}void output(){int i;printf("\n各节点信息:\nxi:");for(i=0;i<N;i++)printf("\t%g",x[i]);printf("\nyi:");for(i=0;i<N;i++)printf("\t%g",y[i]);}结果牛顿插值多项式#include <math.h>#include <stdio.h>int a;#define M 4double x[M+1]={0.4,0.55,0.65,0.8,0.9},y[M+1]={0.41075,0.57815,0.69675,0.88811,1.02652}; void main(){double x;printf("输入x=");scanf("%lf",&x);printf("次数：");scanf("%d",&a);double N(double xx,int a);void output();output();printf("\n%d次牛顿插值多项式N(%g)=%g\n",a,x,N(x,a));}double N(double xx,int a){double s=y[0],d=1;int i,j;double df[M+1][M+1];for(i=0;i<=M;i++)df[i][0]=y[i];for(j=1;j<=a;j++)for(i=j;i<=a;i++)df[i][j]=(df[i][j-1]-df[i-1][j-1])/(x[i]-x[i-j]);printf("\nx\tf(x)\t");for(j=1;j<=a;j++) printf("%5d阶",j);for(i=0;i<=a;i++){printf("\n%g\t%g",x[i],y[i]);for(j=1;j<=i;j++)printf("\t%7.5g",df[i][j]);}for(i=1;i<=a;i++){d*=(xx-x[i-1]);s+=df[i][i]*d;}return s;}void output(){int i;printf("\n各节点信息:\nxi:");for(i=0;i<=M;i++)printf("\t%7g",x[i]);printf("\nyi:");for(i=0;i<=M;i++)printf("\t%7g",y[i]);}结果复合梯形公式#include<stdio.h>#include<math.h>#define f(x) 1/(x*x+1)#define Pi 3.1415926void main(){double a=0,b=1;double T,h,x;int n,i;printf("please input n:");scanf("%d",&n);h=(b-a)/n;x=a;T=0;for(i=1;i<n;i++){x+=h;T+=f(x);}T=(f(a)+2*T+f(b))*h/2;printf("T(%d)=%g\n",n,T);printf("The exact value is %g\n",Pi/4);}复合辛普森公式#include<stdio.h>#include<math.h>#define f(x) 1/(1+x*x)#define Pi 3.1415926void main(){double a=0,b=1;double S,h,x;int n,i;printf("please input Even n:");scanf("%d",&n);h=(b-a)/n;x=a; S=0;for(i=1;i<n;i++){x+=h;if(i%2==0) S+=2*f(x);else S+=4*f(x);}S=(f(a)+S+f(b))*h/3;printf("S(%d)=%g\n",n,S);printf("The exact value is %g\n",Pi/4);}龙贝格公式加速#include<stdio.h>#include<math.h>#define f(x) sin(x)/(1+x)#define M 3void main(){double a=0,b=1;double Long(double a,double b);printf("近似值I=%g\n",Long(a,b));}double Long(double a,double b){int n=1,i=1,j=1;double T[M+1][M+1],h,x,sum;h=b-a;T[0][0]=(f(a)+f(b))/2;for(j=1;j<=3;j++){x=a;h/=2;n*=2;sum=0;for(i=1;i<=n;i+=2){x=a+i*h;sum+=f(x);}T[j][0]=T[j-1][0]/2+h*sum;}for(i=1;i<=M;i++)for(j=1;j<=i;j++){T[i][j]=(pow(4,j)*T[i][j-1]-T[i-1][j-1])/(pow(4,j)-1);}printf("k\tT0\tT1\tT2\tT3\n");for(i=0;i<=M;i++){printf("%d",i);for(j=0;j<=i;j++)printf(" %g",T[i][j]);printf("\n");}return T[M][M];}。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算方法上机报告姓名：学号：班级：上课班级：说明：本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。
1. 对以下和式计算：∑∞⎪⎭⎫ ⎝⎛+-+-+-+=0681581482184161n n n n S n，要求：① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算；（1） 算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 142111416818485861681n n n a n n n n n ε⎛⎫=---<< ⎪+++++⎝⎭； 2、为了保证计算结果的准确性，写程序时，从后向前计算； 3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)（2）算法结构1. ;0=s⎪⎭⎫⎝⎛+-+-+-+=681581482184161n n n n t n； 2. for 0,1,2,,n i =⋅⋅⋅ if 10m t -≤end;3. for ,1,2,,0n i i i =--⋅⋅⋅;s s t =+（3）Matlab源程序clear; %清除工作空间变量clc; %清除命令窗口命令m=input('请输入有效数字的位数m='); %输入有效数字的位数s=0;for n=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));if t<=10^(-m) %判断通项与精度的关系break;endend;fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值for i=n-1:-1:0t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));s=s+t; %求和运算ends=vpa(s,m) %控制s的精度（4）结果与分析当保留11位有效数字时，需要将n值加到n=7，s =3.1415926536；当保留30位有效数字时，需要将n值加到n=22,s =3.14159265358979323846264338328。
通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。
2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:米)如下表所示：①请用合适的曲线拟合所测数据点；②预测所需光缆长度的近似值，作出铺设河底光缆的曲线图；（1）算法思想如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。
分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。
在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。
海底光缆线的长度预测模型如下所示，光缆从A点铺h。
至B点，在某点处的深度为i海底光缆线的长度预测模型计算光缆长度时，用如下公式：20()L f x ds =⎰20'20()1()f x f x dx =+⎰ 191'20()1()k kk f x f x dx+==+∑⎰22()()x y =∆+∆（2）算法结构1. For n i ,,2,1,0⋅⋅⋅=1.1 i i M y ⇒2. For 2,1=k2.1 For k n n i ,,1, -=2.1.1 i k i i i i M x x M M ⇒----)/()(13. 101h x x ⇒-4. For 1-,,2,1n i =4.1 11++⇒-i i i h x x4.2 b a c c h h h i i i i i i ⇒⇒-⇒+++2;1;)/(114.3 i i d M ⇒+165. 0000;;c M d M d n n ⇒⇒⇒λn n n b a b ⇒⇒⇒2;;20μ6. 1111,γμ⇒⇒d b7. 获取M 的矩阵元素个数，存入m 8. For m k ,,3,2 =8.1 k k k l a ⇒-1/μ 8.2 k k k k c l b μ⇒⋅-1- 8.3 k k k k l d γγ⇒⋅-1- 9. m m m M ⇒μγ/10. For 1,,2,1 --=m m k10.1 k k k k k M M c ⇒⋅-+μγ/)(1 11. 获取x 的元素个数存入s 12. k ⇒113. For 1,,2,1-=s i13.1 if i x x ≤~then k i ⇒；break else k i ⇒+114. xx x x x x h x x k k k k ˆ~;~;11⇒-⇒-⇒--- y h x h M y x h M y x M x M k k k k k k ~/]ˆ)6()6(6ˆ6[2211331⇒-+-++---（3）Matlab 源程序clear; clc;x=0:1:20; %产生从0到20含21个等分点的数组X=0:0.2:20;y=[9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7. 95,8.86,9.81,10.80,10.93]; %等分点位置的深度数据n=length(x); %等分点的数目N=length(X);%% 求三次样条插值函数s(x)M=y;for k=2:3; %计算二阶差商并存放在M中for i=n:-1:k;M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));endendh(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量d for i=2:n-1;h(i)=x(i+1)-x(i);c(i)=h(i)/(h(i)+h(i-1));a(i)=1-c(i);b(i)=2;d(i)=6*M(i+1);endM(1)=0; %选择自然边界条件M(n)=0;b(1)=2;b(n)=2;c(1)=0;a(n)=0;d(1)=0;d(n)=0;u(1)=b(1); %对三对角阵进行LU分解y1(1)=d(1);for k=2:n;l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*c(k-1);y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n); %追赶法求解样条参数M(i) for k=n-1:-1:1;M(k)=(y1(k)-c(k)*M(k+1))/u(k);ends=zeros(1,N);for m=1:N;k=1;for i=2:n-1if X(m)<=x(i);k=i-1;break;elsek=i;endendH=x(k+1)-x(k); %在各区间用三次样条插值函数计算X点处的值x1=x(k+1)-X(m);x2=X(m)-x(k);s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-(M(k)*(H^2)/6))*x1+(y(k+1)-(M(k+1)*(H^2)/ 6))*x2)/H;end%% 计算所需光缆长度L=0; %计算所需光缆长度for i=2:NL=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);enddisp('所需光缆长度为L=');disp(L);figureplot(x,y,'*',X,s,'-') %绘制铺设河底光缆的曲线图xlabel('位置','fontsize',16); %标注坐标轴含义ylabel('深度/m','fontsize',16);title('铺设河底光缆的曲线图','fontsize',16);grid;（4）结果与分析铺设海底光缆的曲线图如下图所示：仿真结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图，计算出所需光缆的长度为L=26.4844m。
3. 假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律；试计算这一天的平均气温，并试估计误差。
时刻0 1 2 3 4 5 6 7 8 9 10 11 12平均气温15 14 14 14 14 15 16 18 20 20 23 25 28时刻13 14 15 16 17 18 19 20 21 22 23 24平均气温31 34 31 29 27 25 24 22 20 18 17 16（1）算法思想在本题中，数据点的数目较多。
当数据点的数目很多时，用“多项式插值”方法做数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点是误差很大。
用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数Mi 的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。
另一方面，在有的实际问题中，用插值方法并不合适。
当数据点的数目很大时，要求)(x p 通过所有数据点，可能会失去原数据所表示的规律。
如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。
在这种情况下，不用插值标准而用其他近似标准更加合理。
通常情况下，是选取i α使2E 最小，这就是最小二乘近似问题。
在本题中，采用“最小二乘法”找出这一天的气温变化的规律，使用二次函数、三次函数、四次函数以及指数型函数2)(c t b ae C --=，计算相应的系数，估算误差，并作图比较各种函数之间的区别。
（2）算法结构本算法用正交化方法求数据的最小二乘近似。
假定数据以用来生成了G ，并将y 作为其最后一列（第1+n 列）存放。
结果在a 中，ρ是误差22E 。
I 、使用二次函数、三次函数、四次函数拟合时1. 将“时刻值”存入i x ，数据点的个数存入m2. 输入拟合多项式函数)(x p 的最高项次数1-=n h ，则拟合多项式函数为)()()()(2211x g x g x g x p n n ααα+⋅⋅⋅++= , 根据给定数据点确定G For1,,2,1,0-⋅⋅⋅=n jFor m i ,,2,1⋅⋅⋅=2.1 1,+⇒j i ji g x 2.2 1,+⇒n i i g y3. For n k ,,2,1⋅⋅⋅=3.1 [形成矩阵k Q ]3.1.1 σ⇒-∑=mk i ikkk g g 2/12))(sgn( 3.1.2 k kk g ωσ⇒-3.1.3 For m k k j ,,2,1⋅⋅⋅++=3.1.3.1 j jk g ω⇒ 3.1.4 βσω⇒k 3.2 [变换1-k G 到k G ] 3.2.1 kk g ⇒σFor 1,,,2,1+⋅⋅⋅++=n n k k j 3.2.2 t g mk i ij i ⇒∑=βω/)(3.2.3 For m k k i ,,1,⋅⋅⋅+=3.2.3.1 ij i ij g t g ⇒+ω4. [解三角方程1h Ra =] 4.1nnn n n a g g ⇒+/1.4.2 For 1,,2,1⋅⋅⋅--=n n i 4.2.1iii ni j j ijn i a g x gg ⇒-∑+=+/][11.5. [计算误差22E ]ρ⇒∑+=+mn i n i g121,II 、使用指数函数拟合时现将指数函数进行变形：将y C =,x t =代入2)(c t b aeC --=得:2)(c x b aey --=对上式左右取对数得: 222ln ln bx bcx bc a y -+-=令b bc bc a y z -==-==21202ln ln ααα，，，则可得多项式: 2210x x z ααα++=（3）Matlab 源程序clear; %清除工作空间变量 clc; %清除命令窗口命令 x=0:24; %将时刻值存入数组 y=[15,14,14,14,14,15,16,18,20,20,23,25,28,31,34,31,29,27,25,24,22,20,18,17,16]; [~,m]=size(x); %将数据点的个数存入m T=sum(y(1:m))/m;fprintf('一天的平均气温为T=%f\n',T); %求一天的平均气温 %% 二次、三次、四次函数的最小二乘近似h=input('请输入拟合多项式的最高项次数='); %根据给定数据点生成矩阵G n=h+1; G=[];for j=0:(n-1)g=x.^j; %g(x)按列排列 G=vertcat(G,g); %g 垂直连接G endG=G'; %转置得到矩阵Gfor i=1:m %将数据y 作为G 的最后一列（n+1列）G(i,n+1)=y(i);endG;for k=1:n %形成矩阵Q(k) if G(k,k)>0;sgn=1;elseif G(k,k)==0;sgn=0;else sgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-sgm;for j=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm; %变换Gk-1到Gk for j=k+1:n+1t=sum(w(k:m)*G(k:m,j))/bt;for i=k:m;G(i,j)=G(i,j)+t*w(i);endendendA (n)=G(n,n+1)/G(n,n); %解三角方程求系数Afor i=n-1:-1:1A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);ende=sum(G(n+1:m,n+1).^2); %计算误差efprintf('%d次函数的系数是：',h); %输出系数a及误差e disp(A);fprintf('使用%d次函数拟合的误差是%f：',h,e);t=0:0.05:24;A=fliplr(A); %将系数数组左右翻转Y=poly2sym(A); %将系数数组转化为多项式subs(Y,'x',t);Y=double(ans);figure(1)plot(x,y,'k*',t,Y,'r-'); %绘制拟合多项式函数图形xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');title([num2str(n-1),'次函数的最小二乘曲线']);grid;%% 指数函数的最小二乘近似yy=log(y);n=3;G=[];GG=[];for j=0:(n-1)g=x.^j; %g(x)按列排列G=vertcat(G,g); %g垂直连接Ggg=t.^j; %g(x)按列排列GG=vertcat(GG,gg); %g垂直连接GendG=G'; %转置得到矩阵Gfor i=1:m %将数据y作为G的最后一列（n+1列）G(i,n+1)=yy(i);endG;for k=1:n %形成矩阵Q(k)if G(k,k)>0;sgn=1;elseif G(k,k)==0;sgn=0;else sgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-sgm;for j=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm; %变换Gk-1到Gkfor j=k+1:n+1t=sum(w(k:m)*G(k:m,j))/bt;for i=k:m;G(i,j)=G(i,j)+t*w(i);endendendA(n)=G(n,n+1)/G(n,n); %解三角方程求系数A for i=n-1:-1:1A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);endb=-A(3);c=A(2)/(2*b);a=exp(A(1)+b*(c^2));G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1).^2));e=sum(G(n+1:m,n+1).^2); %计算误差efprintf('\n指数函数的系数是：a=%f,b=%f,c=%f',a,b,c); %输出系数及误差e fprintf('\n使用指数函数拟合的误差是：%f',e);t=0:0.05:24;YY=a.*exp(-b.*(t-c).^2);figure(2)plot(x,y,'k*',t,YY,'r-'); %绘制拟合指数函数图形xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');title(['指数函数的最小二乘曲线']);grid;（4）结果与分析a、二次函数：一天的平均气温为：21.20002次函数的系数: 8.3063 2.6064 -0.0938使用2次函数拟合的误差是：280.339547二次函数的最小二乘曲线如下图所示：b、三次函数：一天的平均气温为：21.20003次函数的系数: 13.3880 -0.2273 0.2075 -0.0084 使用3次函数拟合的误差是：131.061822三次函数的最小二乘曲线如下图所示：c、四次函数：一天的平均气温为：21.20004次函数的系数: 16.7939 -3.7050 0.8909 -0.0532 0.0009 使用4次函数拟合的误差是：59.04118四次函数的最小二乘曲线如下图所示：d、指数函数：一天的平均气温为：21.2000指数函数的系数是：a=26.160286，b=0.004442，c=14.081900使用指数函数拟合的误差是：57.034644指数函数的最小二乘曲线如下图所示：通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果最准确。