计算方法B上机报告
计算方法上机实验报告
. / 《计算方法》上机实验报告班级：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.7（P74）：试利用Jacobi 迭代公式求解方程组要求数值解为方程组的精确解. 二、解题思路：首先将方程组中的系数矩阵A 分解成三部分，即：U D L A ++=，D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。
上机实验报告(精选11篇)
上机实验报告篇1用户名se××××学号姓名学院①实验名称：②实验目的：③算法描述（可用文字描述，也可用流程图）：④源代码：（.c的文件）⑤用户屏幕（即程序运行时出现在机器上的画面）：2．对c文件的要求：程序应具有以下特点：a可读性：有注释。
b交互性：有输入提示。
c结构化程序设计风格：分层缩进、隔行书写。
3.上交时间：12月26日下午1点-6点，工程设计中心三楼教学组。
请注意：过时不候哟！四、实验报告内容0．顺序表的插入。
1.顺序表的删除。
2．带头结点的单链表的\'插入。
3.带头结点的单链表的删除。
注意：1.每个人只需在实验报告中完成上述4个项目中的一个，具体安排为：将自己的序号对4求余，得到的数即为应完成的项目的序号。
例如：序号为85的同学，85%4=1，即在实验报告中应完成顺序表的删除。
2.实验报告中的源代码应是通过编译链接即可运行的。
3.提交到个人空间中的内容应是上机实验中的全部内容。
上机实验报告篇2一、《软件技术基础》上机实验内容1．顺序表的建立、插入、删除。
2．带头结点的单链表的建立（用尾插法）、插入、删除。
二、提交到个人10m硬盘空间的内容及截止时间1．分别建立二个文件夹，取名为顺序表和单链表。
2．在这二个文件夹中，分别存放上述二个实验的相关文件。
每个文件夹中应有三个文件(.c文件、.obj文件和.exe文件)。
3. 截止时间：12月28日（18周周日）晚上关机时为止，届时服务器将关闭。
三、实验报告要求及上交时间（用a4纸打印）1．格式：《计算机软件技术基础》上机实验报告用户名se××××学号姓名学院①实验名称：②实验目的：③算法描述（可用文字描述，也可用流程图）：④源代码：（.c的文件）⑤用户屏幕（即程序运行时出现在机器上的画面）：2．对c文件的要求：程序应具有以下特点：a 可读性：有注释。
b 交互性：有输入提示。
计算方法上机报告_董晓壮
计算方法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 −后形成的矩阵。
数值计算方法上机实验报告
数值计算方法上机实验报告
一、实验目的
本次实验的主要目的是熟悉和掌握数值计算方法，学习梯度下降法的
原理和实际应用，熟悉Python语言的编程基础知识，掌握Python语言的
基本语法。
二、设计思路
本次实验主要使用的python语言，利用python下的numpy，matplotlib这两个工具，来实现数值计算和可视化的任务。
1. 首先了解numpy的基本使用方法，学习numpy的矩阵操作，以及numpy提供的常见算法，如矩阵分解、特征值分解等。
2. 在了解numpy的基本操作后，可以学习matplotlib库中的可视化
技术，掌握如何将生成的数据以图表的形式展示出来。
3. 接下来就是要学习梯度下降法，首先了解梯度下降法的主要原理，以及具体的实际应用，用python实现梯度下降法给出的算法框架，最终
可以达到所期望的优化结果。
三、实验步骤
1. 熟悉Python语言的基本语法。
首先是熟悉Python语言的基本语法，学习如何使用Python实现变量
定义，控制语句，函数定义，类使用，以及面向对象编程的基本概念。
2. 学习numpy库的使用方法。
其次是学习numpy库的使用方法，学习如何使用numpy库构建矩阵，学习numpy库的向量，矩阵操作，以及numpy库提供的常见算法，如矩阵分解，特征值分解等。
3. 学习matplotlib库的使用方法。
计算方法B上机报告
西安交通大学计算方法B上机报告班级：XXXXXXXX姓名：XXX学号：XXXXXXXX2015/12/13 Sunday目录1题目一 (1)1.数值计算 (1)2.实现思想 (1)3.源程序 (1)4.计算结果 (2)5.分析总结 (2)2题目二 (3)1.数据近似 (3)2.实现思想 (3)3.源程序 (5)4.计算结果 (6)5.分析总结 (7)3题目三 (8)1.数据拟合 (8)2.实现思想 (8)3.源程序 (9)4.计算结果 (12)5.分析总结 (14)4题目四 (15)1.非线性方程求解 (15)2.实现思想 (15)3.源程序 (16)4.计算结果 (17)5.分析总结 (17)5题目五 (18)1.线性方程组求解。
(18)2.文件格式说明 (18)3.实现思想 (19)4.源程序 (20)5.计算结果 (23)6.分析总结 (25)7.实际问题 (25)1 题目一1. 数值计算计算以下和式：0142118184858616nn S n n n n ∞=⎛⎫=--- ⎪++++⎝⎭∑，要求：(1)若保留11个有效数字，给出计算结果，并评价计算的算法；(2)若要保留30个有效数字，则又将如何进行计算。
2. 实现思想对于(1)可用迭代的方法进行处理。
通过do-while 循环使和式从0开始计算直到结果满足有效数字即可。
在循环中通过比较本次和上一次结果之差的绝对值与相应有效数字的大小的1/2，即可构成循环中止条件。
对于(2)可用同样的方法进行计算，然而由于所保留的有效数字较大，此时若不对上述算法进行改善的话，由于每一步所得计算结果都是其真实值，此时可能会存在有效数字丢失的情况，从而影响计算的准确性。
因此我们需要对程序中的计算精度进行控制。
在Maltlab 中可利用digits()函数对运算精度进行规定，为防止出现有效位的损失，在实际程序中需将精度提高几位。
而对迭代中每一步所得到的结果则利用vpa()函数进行限定，这样一来我们对每一个运算都控制了精度，而不只是单纯控制了结果的精度3. 源程序使用Matlab 所编写程序如下所示：(1)clear;clc;sd = 11; %所要保留有效数字 n = 0; %迭代次数S1 = 0; %当前n 时计算结果 S2 = 0; %n-1时计算结果 while 1eq = (4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6))/16^n; S1 = S1 + eq;if abs(S1-S2)<0.5*10^(-(sd-1)) %迭代是否终止判断条件 break endS2 = S1; n = n+1; endS=vpa(S1,sd)(2)clear;clc;sd = 30; %所要保留有效数字digits(sd+2);%控制精度n = 0; %迭代次数S1 = vpa(0); %当前n时计算结果S2 = vpa(0); %n-1时计算结果while 1eq = vpa((4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6))/16^n);S1 = vpa(S1 + eq);if abs(S1-S2)<vpa(0.5*10^(-(sd-1))) %迭代是否终止判断条件breakendS2 = S1;n = n+1;endS=vpa(S1,sd)4.计算结果(1)S =3.1415926536(2)S =3.141592653589793238435711679225.分析总结从计算结果可以看出与准确值3.1415926535897932384357116792180407…相比(1)和(2)所得结果分别是准确值保留11位和30位有效数字的结果，即利用本程序所得结果准确。
计算方法上上机实习报告
计算方法上上机实习报告在本次计算方法的上机实习中，我深入体验了数值计算的魅力和挑战，通过实际操作和实践，对计算方法有了更深刻的理解和认识。
实习的目的在于将课堂上学到的理论知识运用到实际的计算中，熟悉各种数值算法的实现过程，提高编程能力和解决实际问题的能力。
我们使用了具体编程语言和软件名称进行编程和计算。
在实习过程中，我首先接触到的是数值逼近的相关内容。
通过多项式插值和曲线拟合的练习，我明白了如何用简单的函数去近似复杂的曲线。
例如，拉格朗日插值法和牛顿插值法让我能够根据给定的离散数据点构建出一个连续的函数，从而对未知点进行预测。
在实际操作中，我需要仔细处理数据的输入和输出，以及算法中的细节，如边界条件和误差控制。
数值积分是另一个重要的部分。
通过梯形公式和辛普森公式，我学会了如何对给定的函数进行数值积分。
在编程实现时，要合理地选择积分区间和步长，以达到所需的精度。
同时，我也了解到了数值积分方法的误差来源和误差估计方法，这对于评估计算结果的可靠性非常重要。
线性方程组的求解是计算方法中的核心内容之一。
我分别使用了高斯消元法和迭代法（如雅克比迭代法和高斯赛德尔迭代法）来求解线性方程组。
在实际编程中，我深刻体会到了算法的效率和稳定性的重要性。
对于大规模的线性方程组，选择合适的算法可以大大提高计算速度和精度。
在非线性方程求根方面，我运用了二分法、牛顿法和割线法等方法。
这些方法各有特点，二分法简单但收敛速度较慢，牛顿法收敛速度快但需要计算导数。
在实际应用中，需要根据方程的特点和求解的要求选择合适的方法。
在实习中，我也遇到了不少问题和挑战。
首先是编程中的错误，如语法错误、逻辑错误等，这需要我耐心地调试和修改代码。
其次，对于一些复杂的算法，理解其原理和实现细节并不容易，需要反复查阅资料和思考。
还有就是数值计算中的误差问题，有时候由于误差的积累，导致计算结果与预期相差较大，需要通过调整算法参数或者采用更精确的算法来解决。
东南大学计算方法上机报告实验报告完整版
实习题11. 用两种不同的顺序计算644834.11000012≈∑=-n n，试分析其误差的变化解：从n=1开始累加，n 逐步增大，直到n=10000；从n=10000开始累加，n 逐步减小，直至1。
算法1的C 语言程序如下： #include<stdio.h> #include<math.h> void main() { float n=0.0; int i; for(i=1;i<=10000;i++) { n=n+1.0/(i*i); } printf("%-100f",n); printf("\n"); float m=0.0; int j; for(j=10000;j>=1;j--) { m=m+1.0/(j*j); } printf("%-7f",m); printf("\n"); }运行后结果如下：结论： 4.设∑=-=Nj N j S 2211，已知其精确值为)11123(21+--N N 。
1）编制按从大到小的顺序计算N S 的程序； 2）编制按从小到大的顺序计算N S 的程序；3）按2种顺序分别计算30000100001000,,S S S ，并指出有效位数。
解：1）从大到小的C语言算法如下：#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void main(){float n=0.0;int i;int N;cout<<"Please input N"<<endl;cin>>N;for(i=N;i>1;i--){n=n+1.0/(i*i-1);N=N-1;}printf("%-100f",n);printf("\n");}执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：2）从小到大的C语言算法如下：#include<stdio.h>#include<math.h>#include<iostream>using namespace std;void main(){float n=0.0;int i;int N;cout<<"Please input N"<<endl;cin>>N;for(i=2;i<=N;i++){n=n+1.0/(i*i-1);}printf("%-100f",n);printf("\n");}执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：结论：通过比较可知：N 的值较小时两种算法的运算结果相差不大，但随着N 的逐渐增大，两种算法的运行结果相差越来越大。
计算方法与计算 实验一误差分析
% 输出的量--每次迭代次数k和迭代值xk,
%
--每次迭代的绝对误差juecha和相对误差xiangcha，
误差分析
误差问题是数值分析的基础，又是数值分析中一个困难的课题。在实际计算 中，如果选用了不同的算法，由于舍入误差的影响，将会得到截然不同的结果。 因此，选取算法时注重分析舍入误差的影响，在实际计算中是十分重要的。同时， 由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差，因此算法 的好坏会影响到数值结果的精度。 一、实验目的
因为运行后输出结果为： y 1.370 762 168 154 49， yˆ =1.370 744 664 189
38， R 1.750 396 510 491 47e-005， WU= 1.782 679 830 970 664e-005 104 . 所
以， yˆ 的绝对误差为 10 4 ，故 y
③ 运行后输出计算结果列入表 1–1 和表 1-2 中。
④ 将算法 2 的 MATLAB 调用函数程序的函数分别用 y1=15-2*x^2 和
y1=x-(2*x^2+x-15)/(4*x+1)代替，得到算法 1 和算法 3 的调用函数程序，将其保
存，运行后将三种算法的前 8 个迭代值 x1, x2 ,, x8 列在一起（见表 1-1），进行
的精确解 x* 2.5 比较，观察误差的传播.
算法 1 将已知方程化为同解方程 x 15 2x2 .取初值 x0 2 ，按迭代公式
xk1 15 2xk2
计算方法与实习上机实验报告
计算方法与实习上机实验报告一、引言本文旨在介绍和展示我们在“计算方法”课程中的实习上机实验环节所完成的一些关键任务和所取得的成果。
该实验课程的目标是让我们更深入地理解和应用各种计算方法，并在实际操作中提高我们的编程和问题解决能力。
二、实验内容与目标实验的主要内容是利用各种计算方法解决实际数学问题。
我们被要求使用编程语言（如Python或Java）来实现和解决这些问题。
这些问题包括使用牛顿法求解平方根，使用蒙特卡洛方法计算圆周率，以及使用最优化方法求解函数的最小值等。
实验的目标不仅是让我们掌握计算方法的基本理论，更是要让我们能够在实际操作中运用这些方法。
我们需要在实习过程中，通过与同伴们合作，共同解决问题，提高我们的团队合作能力和问题解决能力。
三、实验过程与问题解决策略在实验过程中，我们遇到了许多问题，如编程错误、理解困难和时间压力等。
我们通过相互讨论、查阅资料和寻求教师帮助等方式，成功地解决了这些问题。
例如，在实现牛顿法求解平方根时，我们一开始对导数的计算和理解出现了一些错误。
但我们通过查阅相关资料和讨论，最终理解了导数的正确计算方法，并成功地实现了牛顿法。
四、实验结果与结论通过这次实习上机实验，我们不仅深入理解了计算方法的基本理论，还在实际操作中提高了我们的编程和问题解决能力。
我们的成果包括编写出了能有效求解平方根、计算圆周率和求解函数最小值的程序。
这次实习上机实验非常成功。
我们的团队不仅在理论学习和实践操作上取得了显著的进步，还在团队合作和问题解决方面积累了宝贵的经验。
这次实验使我们对计算方法有了更深的理解和认识，也提高了我们的编程技能和解决问题的能力。
五、反思与展望回顾这次实验，我们意识到在实验过程中，我们需要更好地管理我们的时间和压力。
在解决问题时，我们需要更有效地利用我们的知识和资源。
在未来，我们希望能够更加熟练地运用计算方法，并能够更有效地解决问题。
我们也希望能够将所学的计算方法应用到更广泛的领域中，如数据分析、科学研究和工业生产等。
实用数值计算方法上机实验报告
实用数值计算方法上机实验报告学院：化学工程学院姓名：**专业：工业催化学号： **********1. 问题来源某公司饲养实验用的动物以供出售,已知这些动物的生长对饲料中3种营养成分（蛋白质,矿物质和维生素）特别敏感,每个动物每周至少需要蛋白质60g,矿物质3g,维生素8mg,该公司能买到5种不同的饲料,每种饲料1kg 所含各种营养成分和成本如表1所示,如果每个小动物每周食用饲料不超过52kg,求既满足动物生长需要,又能使总成本最低的饲料配方。
数学模型 设需要饲料A1,A2,A3,A4,A5分别为x1,x2,x3,x4,x5（单位kg ）12345min 0.20.70.40.30.5S x x x x x =++++1234512345123451234512350.3x +2x +x +0.6x +1.8x 600.1x +0.05x +0.02x +0.2x +0.05x 3.0.05x +0.1x +0.02x +0.2x +0.08x 8x +x +x +x +x 52,,,4,0s t x x x x x ≥⎧⎪≥⎪⎪≥⎨⎪≤⎪⎪≥⎩在LINGO 的MODEL 窗口内输入如下模型：Min =0.2*x1+0.7*x2+0.4*x3+0.3*x4+0.5*x5; 0.3*x1+2*x2+x3+0.6*x4+1.8*x5>60;0.1*x1+0.05*x2+0.02*x3+0.2*x4+0.05*x5>3; 0.05*x1+0.1*x2+0.02*x3+0.2*x4+0.08*x5>8; x1+x2+x3+x4+x5<52; end求解输出结果如下：Global optimal solution found.Objective value: 22.40000Infeasibilities: 0.000000Total solver iterations: 3Variable Value Reduced Cost X1 0.000000 0.7000000 X2 12.00000 0.000000 X3 0.000000 0.6166667 X4 30.00000 0.000000 X5 10.00000 0.000000 Row Slack or Surplus Dual Price1 22.40000 -1.0000002 0.000000 -0.58333333 4.100000 0.0000004 0.000000 -4.1666675 0.000000 0.8833333结果分析：因此每周每个动物的配料为饲料A2,A4,A5分别为12kg,30kg,10kg,可使得成本达到最低,最低成本为22.4元。
数值计算方法上机实验报告
数值计算方法上机实验报告实验目的：复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。
利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。
上机练习任务：利用计算机基本C 语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。
掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。
一、各算法的算法原理及计算机程序框图1. 列主元高斯消去法算法原理：高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。
列选住院是当高斯消元到第k 步时，从k 列的kk a 以下（包括kk a ）的各元素中选出绝对值最大的，然后通过行交换将其交换到kk a 的位置上。
交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。
●源程序：#define N 200#include "stdio.h"#include "math.h"FILE *fp1,*fp2;void LZ(){int n,i,j,k=0,l;double d,t,t1;static double x[N],a[N][N];fp1=fopen("a1.txt","r");fp2=fopen("b1.txt","w");fscanf(fp1,"%d",&n);for(i=0;i<n;++i)for(j=0;j<=n;++j){fscanf(fp1,"%lf",&a[i][j]);}{d=a[k][k];l=k;i=k+1;do{if(fabs(a[i][k])>fabs(d)) /*选主元*/{d=a[i][k];l=i;}i++;}while(i<n);if(d==0){printf("\n输入矩阵有误！\n");}else{ /*换行*/if(l!=k){for(j=k;j<=n;j++){t=a[l][j];a[l][j]=a[k][j];a[k][j]=t;}}}for(j=k+1;j<=n;j++) /*正消*/ a[k][j]/=a[k][k];for(i=k+1;i<n;i++)for(j=k+1;j<=n;j++)a[i][j]-=a[i][k]*a[k][j];k++;}while(k<n);if(k!=0){for(i=n-1;i>=0;i--) /*回代*/ {t1=0;for(j=i+1;j<n;j++)t1+=a[i][j]*x[j];x[i]=a[i][n]-t1;}for(i=0;i<n;i++)fprintf(fp2,"\n 方程组的根为x[%d]=%lf",i+1,x[i]); fclose(fp1); fclose(fp2); }main() { LZ(); }● 具体算例及求解结果：用列选主元法求解下列线性方程组⎪⎩⎪⎨⎧=++=++=-+28x x 23x 2232832321321321x x x x x x 输入3 输出结果：方程组的根为x[1]=6.0000001 2 -3 8 方程组的根为x[2]=4.000000 2 1 3 22 方程组的根为x[3]=2.000000 3 2 1 28● 输入变量、输出变量说明：输入变量：ij a 系数矩阵元素，i b 常向量元素 输出变量：12,,n b b b 解向量元素2. 杜里特尔分解法解线性方程● 算法原理：求解线性方程组Ax b =时，当对A 进行杜里特尔分解，则等价于求解LUx b =，这时可归结为利用递推计算相继求解两个三角形（系数矩阵为三角矩阵）方程组，用顺代，由Ly b =求出y ，再利用回带，由Ux y =求出x 。
《数值分析》上机实验报告
数值分析上机实验报告《数值分析》上机实验报告1.用Newton 法求方程 X 7-X 4+14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。
1.1 理论依据：设函数在有限区间[a ，b]上二阶导数存在，且满足条件{}αϕ上的惟一解在区间平方收敛于方程所生的迭代序列迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)(')()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20)()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f ab c f x f b a x f b f x f k k k k k k ==-==∈≤-≠>+令)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3225333647>⋅''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f故以1.9为起点⎪⎩⎪⎨⎧='-=+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。
当前后两个的差<=ε时，就认为求出了近似的根。
本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。
1.2 C语言程序原代码：#include<stdio.h>#include<math.h>main(){double x2,f,f1;double x1=1.9; //取初值为1.9do{x2=x1;f=pow(x2,7)-28*pow(x2,4)+14;f1=7*pow(x2,6)-4*28*pow(x2,3);x1=x2-f/f1;}while(fabs(x1-x2)>=0.00001||x1<0.1); //限制循环次数printf("计算结果：x=%f\n",x1);}1.3 运行结果：1.4 MATLAB上机程序function y=Newton(f,df,x0,eps,M)d=0;for k=1:Mif feval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))<=epsd=1;breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey= '奇异'endfunction y=df(x)y=7*x^6-28*4*x^3;Endfunction y=f(x)y=x^7-28*x^4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 问题讨论：1.使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围内得到比较理想的计算结果。
上机报告总结模板下载(通用3篇)
上机报告总结模板下载第1篇一、实验原理（1）电子的自旋轨道磁矩与自旋磁矩 l原子中的电子由于轨道运动，具有轨道磁矩，其数值为：l号表示方向同Pl相反。
在量子力学中PePl2me，负，因而lB1)B2me称为玻尔磁子。
电子除了轨道运动外，其中e还具有自旋运动，因此还具有自旋磁矩，其数值表示为：se Psme。
由于原子核的磁矩可以忽略不计，原子中电子的轨道磁矩和自旋磁矩合成原子的总磁矩：jgej(j1)l(l1)s(s1)Pjg12me，其中g是朗德因子：2j(j1)。
在外磁场中原子磁矩要受到力的作用，其效果是磁矩绕磁场的方向作旋进，也就是Pj 绕着磁场方向作旋进，引入回磁比同时原子角动量Pj和原子总磁矩Pjm ,mj,j1,j2,e2me，总磁矩可表示成jPj。
j取向是量子化的。
Pj在外磁场方向上的投影为：其中m称为磁量子数，相应磁矩在外磁场方向上j。
的投影为： jmmgB ;mj,j1,j2,（2）电子顺磁共振 j。
如果在原子所在的稳定磁场区又叠加一个与之垂直的交变磁场，且角频率满足条件gB B，即EB，刚好满足原子在稳定外磁场中的邻近二能级差时，二邻近能级之间就有共振跃迁，我们称之为电子顺磁共振。
P当原子结合成分子或固体时，由于电子轨道运动的角动量常是猝灭的，即j近似为零，所以分子和固体中的磁矩主要是电子自旋磁矩的贡献。
根据泡利原理，一个电子轨道最多只能容纳两个自旋相反的电子，若电子轨道都被电子成对地填满了，它们的自旋磁矩相互抵消，便没有固有磁矩。
通常所见的化合物大多数属于这种情况，因而电子顺磁共振只能研究具有未成对电子的特殊化合物。
（3）弛豫时间实验样品是含有大量具有不成对电子自旋所组成的系统，虽然各个粒子都具有磁矩，但是在热运动的扰动下，取向是混乱的，对外的合磁矩为零。
当自旋系统处在恒定的外磁场H0中时，系统内各质点的磁矩便以不同的角度取向磁场H0的方向，并绕着外场方向进动，从而形成一个与外磁场方向一致的宏观磁矩M。
计算方法大作业——三次样条插值
计算方法上机报告
此完成所有数据的输入。继续按 Enter 键会出现提示“选择封闭方程组的边界条件: 第 一类边界条件输入 1，第二类边界条件输入 2，第三类边界条件输入 3。 ”根据已知情况 选择相应的边界条件，若为自然三次样条插值，则选 1，并将插值区间两端点的二阶导 数值设置为 0。输入完成之后按 Enter 开始求解，程序运行结束后命令窗口会显示要求 的三次样条插值函数，同时会出现该插值函数以及插值节点的图像，便于直接观察。 2.3 算例及计算结果 (1) 《数值分析》课本第 137 页的例题 4.6.1，已知函数 y=f(x)的数值如下表，求它 的自然三次样条插值函数。 xi yi -3 7 -1 11 0 26 3 56 4 29
(2) 给定函数 f ( x)
3 x 1 1 x 0 0 x3 3 x 4
1 (1 x 1) 。取等距节点，构造牛顿插值多项式 N5(x) 1 25x 2 和 N10(x)及三次样条插值函数 S10(x)。分别将三种插值多项式与 f(x)的曲线画在同一个
N10 x
22757 10 5444 8 20216 6 17147 4 3725 2 x x x x x 1 103 11 53 139 221
将牛顿插值多项式 N5(x)和 N10(x)及三次样条插值函数 S10(x)分别与 f(x)的曲线画在 同一个坐标系上进行比较，如图 12。可以看出三次样条函数与原函数符合的非常好， 对于低次的牛顿插值多项式，与原函数的大致趋势相同，而高次的牛顿插值多项式由 于龙格现象的出现，与原函数之间相差比较大。
S ( xi ) S ( xi ), ( xi ) S ( xi ), S S ( x ) S ( x ), i i i 1, 2, , n 1
电力系统稳态实验报告
电力系统稳态潮流计算上机实验报告一、问题如下图所示的电力系统网络，分别用牛顿拉夫逊法、PQ解耦法、高斯赛德尔法、保留非线性法计算该电力系统的潮流。
发电机的参数如下，*表示任意值负荷参数如下，如上图所示的电力系统，可以看出，节点1、2、3是PQ节点，节点4是PV节点，而将节点5作为平衡节点。
根据问题所需，采用牛顿拉夫逊法、PQ解耦法、高斯赛德尔法、保留非线性法，通过对每次修正量的收敛判据的判断，得出整个电力系统的潮流，并分析这四种方法的收敛速度等等。
算法分析1.牛顿拉夫逊法节点5为平衡节点，不参加整个的迭代过程，节点1、2、3为PQ节点，节点4为PV 节点，计算修正方程中各量，进而得到修正量，判断修正量是否收敛，如果不收敛，迭代继续，如果收敛，算出PQ节点的电压幅值以及电压相角，得出PV节点的无功量以及电压相角，得出平衡节点的输出功率。
潮流方程的直角坐标形式，()()∑∑∈∈++-=ij j ij j ij i ij j ij j ij i i e B f G f f B e G e P()()∑∑∈∈+--=ij j ij j ij i ij j ij j ij i i e B f G e f B e G f Q直角坐标形式的修正方程式，11112n n n m n m -----∆⎡⎤⎡⎤∆⎡⎤⎢⎥⎢⎥∆=-⎢⎥⎢⎥⎢⎥∆⎣⎦⎢⎥⎢⎥∆⎣⎦⎣⎦PHN e Q M L f UR S修正方程式中的各量值的计算，()()][∑∑∈∈++--=∆ij j ij j ij i ij j ij j ij i is i e B f G f f B e G e p P()()][∑∑∈∈+---=∆ij j ij j ij i ij j ij j ij i is i e B f G e f B e G f Q Q)(2222i i is i f e U U +-=∆Jacobi 矩阵的元素计算，()()()ij i ij i i ijij j ij j ii i ii i jj iB e G f j i Q M G f B e B e G f j i e ∈-⎧≠∂∆⎪==⎨++-=∂⎪⎩∑()()()ij i ij i i ijij j ij j ii i ii i jj iG e B f j i Q L G e B f G e B f j i f ∈+⎧≠∂∆⎪==⎨--++=∂⎪⎩∑)()(202i j i j e e U R ijij i =≠⎩⎨⎧-=∂∆∂=)()(202i j i j f f U S ijij i =≠⎩⎨⎧-=∂∆∂=牛顿拉夫逊法潮流计算的流程图如下，2.PQ 解耦法如同牛顿拉夫逊法，快速解耦法的前提是，输电线路的阻抗要比电阻大得多，并且输电线路两端的电压相角相差不大，此时可利用PQ 快速解耦法，来计算整个电力系统网络的潮流。
电力系统潮流计算实验报告
１1． 手算过程已知：节点1：PQ 节点， s(1)= -0.5000-j0.3500 节点2：PV 节点， p(2)=0.4000 v(2)=1.0500 节点3：平衡节点，U(3)=1.0000∠0.0000 网络的连接图：0.0500+j0.2000 1 0.0500+j0.2000231）计算节点导纳矩阵由2000.00500.012j Z 71.418.112j y ;2000.00500.013j Z71.418.113j y ;导纳矩阵中的各元素：42.936.271.418.171.418.1131211j j j y y Y ;71.418.11212j y Y ; 71.418.11313j y Y; 21Y 71.418.11212j y Y ; 71.418.12122j y Y;002323j y Y;31Y 71.418.11313j y Y; 32Y 002323j y Y;71.418.13133j y Y;形成导纳矩阵BY ：71.418.10071.418.10071.418.171.418.171.418.171.418.142.936.2j j j j j j j j j Y B2）计算各PQ、PV 节点功率的不平衡量，及PV 节点电压的不平衡量：取:000.0000.1)0(1)0(1)0(1j jf e U000.0000.1)0(2)0(2)0(2j jf e U节点3是平衡节点，保持000.0000.1333j jf e U为定值。
nj j jij jij ijij jij i ieB fG f fB eG e P1)0()0()0()0()0()0()0(；２nj j jij jij ijij jij i ie B fG e f B eG f Q 1)0()0()0()0()0()0()0(；);(2)0(2)0(2)0(iiif e U)0.142.90.036.2(0.0)0.042.90.136.2(0.1)0(1P)0.171.40.018.1(0.0)0.071.40.118.1(0.1 )0.171.40.018.1(0.0)0.071.40.118.1(0.1 0.0 ；)0.142.90.036.2(0.1)0.042.90.136.2(0.0)0(1Q)0.171.40.018.1(0.1)0.071.40.118.1(0.0 )0.171.40.018.1(0.1)0.071.40.118.1(0.0 0.0 ；)0.171.40.018.1(0.0)0.071.40.118.1(0.1)0(2P)0.171.40.018.1(0.0)0.071.40.118.1(0.1 )0.00.00.00.0(0.0)0.10.00.10.0(0.1 0.0 ；101)(222)0(22)0(22)0(2f e U；于是：;)0()0(iiiP P P ;)0()0(iiiQQ Q);(2)0(2)0(22)0(iiiif e UU5.00.05.0)0(11)0(1P P P ;35.00.035.0)0(11)0(1QQ Q;4.00.04.0)0(22)0(2P P P ;1025.0)01(05.1)(2222)0(22)0(2222)0(2f e UU3）计算雅可比矩阵中各元素雅可比矩阵的各个元素分别为：３ji ij ji ij j i ij j i ij ji ij j i ij e U S f U R e Q L f Q J e P N f P H 22;;; 又： nj j jij jij i jij jij i ieB fG f fB eG e P1)0()0()0()0()0()0()0(； nj j jij jij ijij jij iieB fG e fB eG f Q 1)0()0()0()0()0()0()0(；);(2)0(2)0(2)0(iiif e U)0(1P )0(111)0(111)0(1)0(111)0(111)0(1e Bf G f f B e G e)0(212)0(212)0(1)0(212)0(212)0(1e B fG f f B e G e313313)0(1313313)0(1e Bf G f f B e G e ；)()()0(111)0(111)0(1)0(111)0(111)0(1)0(1e Bf Ge f B e G f Q)()()0(212)0(212)0(1)0(212)0(212)0(1e Bf G e f B e G f)()(313313)0(1313313)0(1e Bf G e f B e G f；)0(2P )0(121)0(121)0(2)0(121)0(121)0(2e Bf G f f B e G e)0(222)0(222)0(2)0(222)0(222)0(2eB fG f fBeG e323323)0(2323323)0(2e Bf G f f B e G e ；)(2)0(22)0(22)0(2f e U42.90.171.40.171.4313)0(212)0(1)0(1)0(11e B e Bf P H ； 36.20.118.10.118.10.136.222313)0(212)0(111)0(1)0(1)0(11 e G e G e G e P N 36.20.118.10.118.1313)0(212)0(1)0(1)0(11 e G e G f Q J４42.90.171.40.171.40.142.922313)0(212)0(111)0(1)0(1)0(11 e B e B e B e Q L 71.40.171.4)0(112)0(2)0(1)0(12 e B f P H ； 18.10.118.1)0(112)0(2)0(1)0(12 e G e P N ； 18.10.118.1)0(112)0(2)0(1)0(12 e G f Q J ；71.40.171.4)0(112)0(2)0(1)0(12 e B e Q L ； 71.40.171.4)0(221)0(1)0(2)0(21 e B f P H ； 11.40.111.4)0(221)0(1)0(2)0(21 e G e P N ； 0)0(12)0(2)0(21 f U R ； 0)0(12)0(2)0(21 e U S ； 71.40.10.00.171.4323)0(121)0(2)0(2)0(22 e B e B f P H ； 18.10.10.00.118.10.118.122323)0(121)0(222)0(2)0(2)0(22 e G e G e G e P N ；02)0(2)0(22)0(2)0(22 f f U R ； 0.20.122)0(2)0(22)0(2)0(22 e e U S ； 得到K=0时的雅可比矩阵：0.200018.171.418.171.471.418.142.936.218.171.436.242.9)0(J4）建立修正方程组：５)0(2)0(2)0(1)0(10.200011.4959.1011.4959.10959.1011.4918.2122.811.4959.1022.8918.210975.04.035.08.0e f e f 解得：04875.001828.00504.00176.0)0(2)0(2)0(1)0(1e f e f 因为 )0()0()1(iiie e e ； )0()0()1(iiif f f ；所以 9782.00218.00.1)0(1)0(1)1(1e e e ； 0158.00158.00)0(1)0(1)1(1f f f ；05125.105125.00.1)0(2)0(2)1(2e e e ；05085.005085.00)0(2)0(2)1(2f f f ；5）运用各节点电压的新值进行下一次迭代：即取: 0158.09782.0)1(1)1(1)1(1j jf e U05085.005125.1)1(2)1(2)1(2j jf e U节点3时平衡节点，保持000.0000.1333j jf e U为定值。
计算方法实验上机报告(完整版)
简单迭代法#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)。
计算方法B上机实习报告计算方法B 上机实习报告1.对以下和式计算：0142118184858616n n S n n n n ∞=⎛⎫=--- ⎪++++⎝⎭∑，要求： (1)若只需保留11个有效数字，该如何进行计算；(2)若要保留30个有效数字，则又将如何进行计算；问题分析：在该题中S 的每一项14211()1681848586n ns n n n n =---++++存在两个相近的数相减的问题，因此为了避免有效数字损失，最好是改变运算顺序，分别将正数和负数相加，然后再将其和相加。
另外，sn 中有多个负数相加，可以按照绝对值递增的顺序求和，以减少舍入误差的影响。
同时，为了避免大数吃小数的问题，本题先计算出保留目标有效数字所需 要的迭代次数，然后采用倒序相加的方法实现。
程序实现：clear;clc;m=input('请输入要保留的有效数字位数：'); s1=0; s2=0; k=0; s=1;%%%%判断多需要的迭代次数 while s>=0.5*10^-(m-1)s=4/(16^k*(8*k+1))-(2/(16^k*(8*k+4))+1/(16^k*(8*k+5))+1/(16^k*(8*k+6))); k=k+1; end%%%%正负数分别按照绝对值递增的顺序倒序相加 for n=(k-1):-1:0 a1=4/(16^n*(8*n+1)); a2=2/(16^n*(8*n+4)); a3=1/(16^n*(8*n+5)); a4=1/(16^n*(8*n+6)); s1=a1+s1; s2=a4+a3+a2+s2; end S=s1-s2; S=vpa(S,m)运算结果：总结心得：在计算求和问题中，应特别注意相近数相减的问题，这样会造成有效数字灾难性的损失。
另外在两个数量级相差较大的数字相加减时，较小数的有效数字会被丧失，因此要按照从小到大的顺序相加。
在上题计算中分别对正负相采用倒序相加，这样就有效的避免了“大数吃小数”的问题。
2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:米)(1)请用合适的曲线拟合所测数据点；(2)预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；问题分析：本题的主要目的是对测量数据进行拟合，同时对拟合曲线进行线积分即可得到河底光缆长度的近似值。
由于数值点较多时，使用拉格朗日差值多项式会出现龙格现象。
为了将所有的数据点都用上，采用分段差插法，本题使用三次样条插值。
算法思想：M记样条函数在每个子区间上是三次多项式，它的二阶导数必是一次多项式。
若用i在i x 处的二阶导数''(x )i s 。
则在区间1[x ,x ]i i - 上''11(x)i i i i ix x x x s M M h h ----=+ 式中 1i i i h x x -=- （1） 对上式进行两次积分得332211111()()()()()6666i i i i i i i i i i i i i i i ix x x x h x x h x x s x M M y M y M h h h h ---------=++-+- （2） 它的一阶导数为22'1111()()()226i i i i i i i i i i i i x x x x y y M M s x M M h h h h --------=-++- （3）'()i s x 满足连续性条件，即''(x )()i i s s x -+= 1,2,n i =-…,1上式和（3）式得1111111636i i i i i i i i i i i i ih h h h y y y y M M M h h +++--+++--++=- （4） 用差商记号，并记11i i i i h h h λ++=+ 1ii i i h h h μ+=+（4）式可以写成111126[,,]i i i i i i i i M M M y x x x μλ-+-+++= 1,2,1i n =-…,方程组可以写成如下形式00111111112222n n n n nn n M d M d M d M d λμλμλμ----⎛⎫⎛⎫⎛⎫⎪⎪ ⎪ ⎪⎪ ⎪ ⎪⎪ ⎪= ⎪⎪ ⎪ ⎪⎪ ⎪⎪⎪ ⎪⎝⎭⎝⎭⎝⎭O OO M M 116[,,]i i i i d y x x x -+= i 1,2,n =…,自然样条插值条件为000,0,0,0n n d d λμ====在估计河底光缆长度时使用第一类线积分2019k kk l ds +====∑⎰⎰⎰长度程序实现：clear;clc;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.80 10.93];d=y;plot(x,y,'k.','markersize',15)hold on%%%计算牛顿二阶差商for k=1:2for i=21:-1:(k+1)d(i)=(d(i)-d(i-1))/(x(i)-x(i-k));endend%%%假定d的边界条件，采用自然三次样条for i=2:20d(i)=6*d(i+1);endd(1)=0;d(21)=0;%%%追赶法求解带状矩阵的m值a=0.5*ones(1,21);b=2*ones(1,21);c=0.5*ones(1,21);a(1)=0;c(21)=0;u=ones(1,21);u(1)=b(1);r=c;yy(1)=d(1);%%%追的过程for k=2:21l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*r(k-1);yy(k)=d(k)-l(k)*yy(k-1);end%%%赶的过程m(21)=yy(21)/u(21);for k=20:-1:1m(k)=(yy(k)-r(k)*m(k+1))/u(k);end%%%利用插值点画出拟合曲线k=1;nn=100;xx=linspace(0,20,nn);l=0;for j=1:nnfor i=2:20if xx(j)<=x(i)k=i;break;elsek=i+1;endendh=1;xbar=x(k)-xx(j);xmao=xx(j)-x(k-1);s(j)=(m(k-1)*xbar^3/6+m(k)*xmao^3/6+(y(k-1)-m(k-1)*h^2/6)*xbar+(y(k)-m(k)*h^2/6)*xmao)/h; sp(j)=-m(k-1)*(x(k)-xx(j))^2/(2*h)+m(k)*(xx(j)-x(k-1))^2/(2*h)+(y(k)-y(k-1))/h-(m(k)-m(k-1))*h/6; l(j+1)=(1+sp(j)^2)^0.5*(20/nn)+l(j);%利用第一类线积分求河底光缆的长度end%%%绘图plot(xx,s,'r-','linewidth',1.5)griddisp(['所需光缆长度为',num2str(l(nn+1)),'米'])运行结果：总结心得：采用三次样条插值对数据进行拟合时，可以有效避免龙格现象。
在本题的计算中采用自然三次样条函数的边界条件。
在解线性方程组时使用了追赶法求解带状矩阵，在求解三对角矩阵时追赶法计算速度快，是一种求解线性方程组的有效手段。
在估计河底光缆长度时使用第一类线积分。
本题计算中间变量非常多，在调试的过程中遇到了一些麻烦，这更加使我认识到在编程的过程中由不得一点儿马虎，在下面问题的处理中我变得更加小心。
3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温问题分析：对一天的气温进行数据拟合，可以考虑使用最小二乘的二次函数、三次函数、四次函数以及指数函数。
问题的难度是求解各种拟合函数的系数。
利用法方程求解最小二乘系数时，方程的解不够稳定，本题利用正交化方法。
程序实现：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=length(x);n=input('请输入函数的次数');plot(x,y,'k.',x,y,'-')grid;hold on;n=n+1;G=zeros(m,n+1);G(:,n+1)=y';c=zeros(1,n);%建立c来存放σq=0;f=0;b=zeros(1,m);%建立b用来存放β%%%形成矩阵Gfor j=1:nfor i=1:mG(i,j)=x(1,i)^(j-1);endend%%%建立矩阵Qkfor k=1:nfor i=k:mc(k)=G(i,k)^2+c(k);endc(k)=-sign(G(k,k))*(c(k)^0.5);w(k)=G(k,k)-c(k);%建立w来存放ωfor j=k+1:mw(j)=G(j,k);endb(k)=c(k)*w(k);%%%变换矩阵Gk-1到GkG(k,k)=c(k);for j=k+1:n+1q=0;for i=k:mq=w(i)*G(i,j)+q;ends=q/b(k);for i=k:mG(i,j)=s*w(i)+G(i,j);endendend%%%求解三角方程Rx=h1a(n)=G(n,n+1)/G(n,n);for i=n-1:(-1):1for j=i+1:nf=G(i,j)*a(j)+f;enda(i)=(G(i,n+1)-f)/G(i,i); %a（i）存放各级系数f=0;enda%%%拟合后的回代过程p=zeros(1,m);for j=1:mfor i=1:np(j)=p(j)+a(i)*x(j)^(i-1);end endplot(x,p,'r*',x,p,'-'); E2=0;%用E2来存放误差 %%%误差求解 for i=n+1:mE2=G(i,n+1)^2+E2; end E2=E2^0.5; disp('误差为'); disp(E2); t=0 for i=1:m t=t+p(i); endt=t/m; %%%平均温度disp(['平均温度为',num2str(t),'℃'])运行结果：二次函数时，系数a1=8.3063,a2=2.6064,a3=-0.0938。