计算方法上机实践

计算方法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 −后形成的矩阵。

计算方法上机实习报告[题目及目的要求]1.用二分法求方程0163=--x x 在[0,5]上根的近似值。

用牛顿迭代法求0133=--x x 在2=x 附近的实根。

2.完成史蒂芬森迭代加速法和割线法的子程序，并利用方程010423=--x x 对比对分法与一般迭代法。

[方法原理说明]1.二分法和牛顿迭代法：二分法是逐次把有根区间分半，舍弃无根区间而保留有根区间的一种逼近根的方法。

在这个过程中有根区间的长度以2的幂次方减少，当有根区间的长度小于给定的精度时，其中点就作为根的近似值。

牛顿迭代法的迭代格式为： 初值0x ()()k k k k x f x f x x '1-=+ （k=0,1,2....）显然，牛顿迭代格式能够迭代下去必须要求()k x f 的导数不能为0.当某个()0'=k x f 或很小时，迭代中断；当()k x f 满足一定条件时，牛顿迭代具有平方收敛速度。

该方法对初值0x 要求较高，若选取不当，则可能发散，若选取的好，则收敛很快。

2.史蒂芬森迭代加速法和割线法迭代法就是通过一个迭代格式进行反复迭代以产生一个序列。

若这个序列收敛于方程的根，就称这个迭代格式收敛。

史蒂芬森迭代加速法的迭代格式为()k k k k k k k x y z x y x x +---=+221()k k x f y = ，()k k y f z = （k=0,1,2....）割线法与牛顿迭代法一样，即在根的某个邻域内，()k x f 有直至二阶的连续导数，且()0'≠k x f ，则在邻域内选取初值10,x x ，迭代均收敛。

割线法的迭代格式为初值10,x x()()()()111--+---=k k k k k k k x x x f x f x f x x （k=2,3....）[计算步骤]1.二分法:1）给定a,b 及精度要求ep ； 2）计算x=（a+b ）/2 及()k x f ；3）若b-a<ep ，则返回主程序，x 作为近似根，否则转4； 4）若()()0<a f x f ，则b x ⇒，否则a x ⇒； 5）转2。

计算方法上机作业——龙格库塔法龙格库塔法（Runge-Kutta method）是一种常用于求解常微分方程（Ordinary Differential Equation，ODE）的数值解法。

它是由德国数学家卡尔·龙格（Carl Runge）和马丁·威尔海姆·库塔（Martin Wilhelm Kutta）在20世纪初提出的。

龙格库塔法的基本思想是通过数值逼近来计算微分方程的近似解。

在讲解龙格库塔法之前，我们先来简单回顾一下ODE的一阶常微分方程的基本形式：y′(y)=y(y,y)其中，y(y,y)是已知函数。

龙格库塔法的核心是使用差分逼近计算函数的斜率。

假设我们要求解的方程为：y′(y)=y(y,y),y(y)=y₀所需计算的点为y₀,y₁,...,yy，对应的函数值为y₀,y₁,...,yy，其中y是步长的个数。

龙格库塔法通过递推关系式来计算估计值，并不断更新当前点的函数值。

接下来以龙格库塔法的经典四阶形式为例进行说明。

该方法的基本方程如下：yy+1=yy+(y₁+2y₂+2y₃+y₄)/6y₁=ℎy(yy,yy)y₂=ℎy(yy+ℎ/2,yy+y₁/2)y₃=ℎy(yy+ℎ/2,yy+y₂/2)y₄=ℎy(yy+ℎ,yy+y₃)其中y表示当前步骤，ℎ表示步长，yy表示当前点的函数值，y₁,y₂,y₃和y₄则表示对应的斜率。

使用龙格库塔法，我们可以通过不断递归计算来求得指定区间（例如[y,y]）上的函数值。

具体步骤如下：1.确定求解区间[y,y]和初始点(y₀,y₀)以及步长ℎ。

2.初始化：设置yy=y₀，yy=y₀。

3.对所有y=0,...,y−1：计算y₁,y₂,y₃和y₄，根据上述递推关系式。

根据递推关系式计算yy+1更新当前点的函数值，即yy+1=y(yy+1)。

更新当前点的y值，即yy+1=yy+ℎ。

4.返回结果：最终求得的函数值。

需要注意的是，选择适当的步长对最终结果的精度和计算效率都有重要影响。

计算方法上机实习题目（四）——龙贝格算法计算椭圆周长一、 题目用龙贝格算法计算椭圆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]，我在最初编写程序时就忽略了这个问题，导致花费很久调试程序。

数值计算方法上机实验报告
一、实验目的
本次实验的主要目的是熟悉和掌握数值计算方法，学习梯度下降法的
原理和实际应用，熟悉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库的使用方法。

东南大学计算方法与实习实验报告学院：电子科学与工程学院学号：06A*****姓名：***指导老师：***实习题14、设S N=Σ(1)编制按从大到小的顺序计算S N的程序；(2)编制按从小到大的顺序计算S N的程序；(3)按两种顺序分别计算S1000，S10000，S30000，并指出有效位数。

解析：从大到小时，将S N分解成S N-1=S N-,在计算时根据想要得到的值取合适的最大的值作为首项；同理从小到大时，将S N=S N-1+ ，则取S2=1/3。

则所得式子即为该算法的通项公式。

(1)从大到小算法的C++程序如下：/*从大到小的算法*/#include<iostream>#include<iomanip>#include<cmath>using namespace std;const int max=34000; //根据第(3)问的问题，我选择了最大数为34000作为初值void main(){int num;char jus;double cor,sub;A: cout<<"请输入你想计算的值"<<'\t';cin>>num;double smax=1.0/2.0*(3.0/2.0-1.0/max-1.0/(max+1)),temps;double S[max];// cout<<"s["<<max<<"]="<<setprecision(20)<<smax<<'\n';for(int n=max;n>num;){temps=smax;S[n]=temps;n--;smax=smax-1.0/((n+1)*(n+1)-1.0);}cor=1.0/2.0*(3.0/2.0-1.0/num-1.0/(num+1.0)); //利用已知精确值公式计算精确值sub=fabs(cor-smax); //double型取误差的绝对值cout<<"用递推公式算出来的s["<<n<<"]="<<setprecision(20)<<smax<<'\n';cout<<"实际精确值为"<<setprecision(20)<<cor<<'\n';cout<<"则误差为"<<setprecision(20)<<sub<<'\n';cout<<"是否继续计算S[N],是请输入Y，否则输入N！"<<endl;cin>>jus;if ((int)jus==89||(int)jus==121) goto A;}(2)从小到大算法的C++程序如下：/*从小到大的算法*/#include<iostream>#include<iomanip>#include<cmath>using namespace std;void main(){int max;A: cout<<"请输入你想计算的数，注意不要小于2"<<'\t';cin>>max;double s2=1.0/3.0,temps,cor,sub;char jus;double S[100000];for(int j=2;j<max;){temps=s2;S[j]=temps;j++;s2+=1.0/(j*j-1.0);}cor=1.0/2.0*(3.0/2.0-1.0/j-1.0/(j+1.0)); //利用已知精确值公式计算精确值sub=fabs(cor-s2); //double型取误差的绝对值cout<<"用递推公式算出来的s["<<j<<"]="<<setprecision(20)<<s2<<'\n';cout<<"实际精确值为"<<setprecision(20)<<cor<<'\n';cout<<"则误差为"<<setprecision(20)<<sub<<'\n';cout<<"是否继续计算S[N],是请输入Y，否则输入N！"<<endl;cin>>jus;if ((int)jus==89||(int)jus==121) goto A;}(3)(注：因为程序中setprecision(20)表示输出数值小数位数20，则程序运行时所得到的有效数字在17位左右)ii.选择从小到大的顺序计算S1000、S10000、S30000的值需要计算的项S1000S10000S30000计算值0.74900049950049996 0.74966672220370571 0.74996666722220728实际精确值0.74900049950049952 0.74990000499950005 0.74996666722220373误差 4.4408920985006262*10-16 5.6621374255882984*10-15 3.5527136788005009*10-15有效数字17 17 17附上部分程序运行图：iii.实验分析通过C++程序进行计算验证采用从大到小或者从小到大的递推公式算法得到的数值基本稳定且误差不大。

实习题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 的逐渐增大，两种算法的运行结果相差越来越大。

数值分析上机实践报告一、实验目的本次实验主要目的是通过上机操作，加深对数值分析算法的理解，并熟悉使用Matlab进行数值计算的基本方法。

在具体实验中，我们将实现三种常见的数值分析算法：二分法、牛顿法和追赶法，分别应用于解决非线性方程、方程组和线性方程组的求解问题。

二、实验原理与方法1.二分法二分法是一种常见的求解非线性方程的数值方法。

根据函数在给定区间端点处的函数值的符号，不断缩小区间的长度，直到满足精度要求。

2.牛顿法牛顿法是求解方程的一种迭代方法，通过构造方程的泰勒展开式进行近似求解。

根据泰勒展式可以得到迭代公式，利用迭代公式不断逼近方程的解。

3.追赶法追赶法是用于求解三对角线性方程组的一种直接求解方法。

通过构造追赶矩阵，采用较为简便的向前追赶和向后追赶的方法进行计算。

本次实验中，我们选择了一组非线性方程、方程组和线性方程组进行求解。

具体的实验步骤如下：1.调用二分法函数，通过输入给定区间的上下界、截止误差和最大迭代次数，得到非线性方程的数值解。

2.调用牛顿法函数，通过输入初始迭代点、截止误差和最大迭代次数，得到方程组的数值解。

3.调用追赶法函数，通过输入追赶矩阵的三个向量与结果向量，得到线性方程组的数值解。

三、实验结果与分析在进行实验过程中，我们分别给定了不同的参数，通过调用相应的函数得到了实验结果。

下面是实验结果的汇总及分析。

1.非线性方程的数值解我们通过使用二分法对非线性方程进行求解，给定了区间的上下界、截止误差和最大迭代次数。

实验结果显示，根据给定的输入，我们得到了方程的数值解。

通过与解析解进行比较，可以发现二分法得到的数值解与解析解的误差在可接受范围内，说明二分法是有效的。

2.方程组的数值解我们通过使用牛顿法对方程组进行求解，给定了初始迭代点、截止误差和最大迭代次数。

实验结果显示，根据给定的输入，我们得到了方程组的数值解。

与解析解进行比较，同样可以发现牛顿法得到的数值解与解析解的误差在可接受范围内，说明牛顿法是有效的。

数值计算方法第一次实习报告一、 实习题目。

1. 分别用高斯-赛德尔迭代法、雅克比迭代法、列主元消去法解下列方程组。

⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++=+++ 1.84711.2671x 0.2568x 0.2471x 0.2368x1.74710.2271x 1.2168x 0.2071x 0.1968x 1.64710.1871x 0.1768x 1.1675x 0.1582x 1.54710.1490x 0.1397x 0.1254x 1.1161x4321432143214321 2. 用迭代法求x 5-x-0.2=0的正根，要求精确到小数点后第五位。

二、算法原理。

1、雅可比迭代法基本原理 将矩阵分解为，其中则式可记为，变形可得，可逆时，有于是得到迭代的过程为式中，，即2、高斯-赛德尔迭代法基本原理赛德尔迭代法是对雅可比迭代法的一种改进，雅可比迭代法是在每一步计算的各个分量时均只用到中的分量。

实际上，在计算时，分量都已经计算出来而没有被直接利用，因此可以考虑以来代替计算。

即121311121232222121-1,212121000, ,0000n n n n nn a a a a a a a a D L U aa a a a a a ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==-=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦=Ax b ()=D L U x b--()=+Dx L U x b+D ()11220nn a a a ≠ ()11=D +x L U x D b --+=+x Bx f ()11D,=D B L U f b --=+()=11=+=1,2,,n i i ij j j ii j i x b a x i n a ≠⎛⎫ ⎪ ⎪ ⎪⎝⎭∑ ()+1k x ()k x ()+1k i x ()()+1+1-1,,k k i i x x ()()+1+1-1,,k k i i x x ()()1-1,,kki x x矩阵形式为，可得，于是赛德尔迭代法的矩阵形式为式中，。