数值计算课程设计实验合并

1、本科生实验报告实验课程 数值计算方法 学院名称 信息科学与技术学院 专业名称 计算机科学与技术 学生姓名 周瑞豪 学生学号 202113030317 指导教师 马桂媛 实验地点 6A502 实验成绩 二 15 年 4 月 二 15 年 5 月填写说明1、 适用于本科生所有的实验报告印制实验报告册除外；2、 专业填写为专业全称，有专业方向的用小括号标明；3、 格式要求： 用A4纸双面打印封面双面打印或在A4大小纸上用蓝黑色水笔书写。 打印排版：正文用宋体小四号，1.5倍行距，页边距采取默认形式上下2.54cm，左右2.54cm，页眉1.5cm，页脚1.75cm。字符间距为默认值缩放100%，间距

2、：标准；页码用小五号字底端居中。 具体要求：题目二号黑体居中；摘要“摘要二字用小二号黑体居中，隔行书写摘要的文字局部，小4号宋体；关键词隔行顶格书写“关键词三字，提炼3-5个关键词，用分号隔开，小4号黑体)； 正文局部采用三级标题；第1章 ××(小二号黑体居中，段前0.5行)1.1 ×××××小三号黑体×××××段前、段后0.5行小四号黑体段前、段后0.5行参考文献黑体小二号居中，段前0.5行，参考文献用五号宋体，参照?参考文献著录规那么GB/T 77142005?。实验一

3、非线性方程求根第1章 问题描述实验目的：掌握非线性方程求根的根本步骤及方法，。实验内容：试分别用二分法、简单迭代法、Newton迭代法、弦截法(割线法、双点弦法)，求 x5-3x3+x-1= 0 在区间 -8,8上的全部实根，误差限为10-6。要求：讨论求解的全过程，对所用算法的局部收敛性，优缺点等作分析及比拟， 第2章 算法思想 算法分析此局部摘抄自课本P23页 设函数有根区间表示为。将有根区间用中点分成两 半，计算函数值。如果，就得到方程组的实根，否那么检查与是否同号，假设同号，那么说明所求的根在的右侧，这时令；否那么，根在的左侧，这时令，这样新的有根区间的长度为之半。对压缩了的有根区间又

4、施以同样的方法如此反复二分下去，即可得出一系列有根区间其中，每个区间都是前一个区间的一半，因此二分k次后的有根区间的长度为可见，如果二分过程无限地下去，这些有根区间最终必收缩于一点，该点显然就是所求的根。取有根区间的中点作为根的近似值，此时的误差假设事先给定的误差要求为，那么只需便可以停止二分计算。 void divide(V a1,V a2) for(;fabs(fabs(a1)-fabs(a2)>=1e-6;) if(f(a1)*f(a1+a2)/2)<0) a2=(a1+a2)/2;else a1=(a1+a2)/2; cout<<endl<<&quo

5、t;经过二分法求得的结果是："<<(a1+a2)/2<<endl<<endl; 2.2 简单迭代法 算法分析此局部摘抄自肯课本P41面Newton下山法是扩大初值范围的修正Newton法。将Newton迭代法的计算结果进行适当的加权平均作为新的改良值，即。从而化简可得简单法迭代公式 void newton(V a1) for(;fabs(fabs(a1)-fabs(a1-fnewton(a1)>=1e-6;) a1-=fnewton(a1); cout<<"经过简单迭代法得的结果是："<<a1<

6、;<endl<<endl; 2.3 Newton迭代法 算法分析(此局部摘抄自书上P37面) 对于非线性方程，假设根的一个近似值在这里是，将在处展成一阶泰勒公式，忽略高次项，有。右端是直线方程，用这个直线方程来近似非线性方程。将非线性方程的根代入，即解得这就是Newton迭代公式。在具体的应用中，写为：，那么这样获得的即为按Newton迭代法求得的近似解。 void newton2(V a1) double t=1.0; for(;fabs(f(a1)<=fabs(f(a1-t*fnewton(a1);t/=2.0) for(;fabs(fabs(a1)-fabs(a1

7、-t*fnewton(a1)>=1e-6;) a1-=t*fnewton(a1); cout<<"使用Newton迭代法求得的结果是："<<a1<<endl<<endl; 算法分析此局部摘抄自课本P44面由于Newton迭代法有时收敛速度较慢，而且有时函数的一阶导数不易求得或较为复杂，因此，改用两个端点都在变动的弦，即用差商替代Newton迭代公式中的导数，从而导出。这就是双点弦截法。双点弦截法详细的计算步骤为：1选定初始值，计算，。2按双点弦法迭代公式计算，并求。3判断：如果给定精度，那么迭代停止；否那么，用和分别代替

8、重复2和3。第3章 测试结果及分析 测试结果，要求按表格形式输出每次迭代的结果表格格式如书上例题，并比拟不同方法的收敛速度、优劣附录：源程序清单#include ""#include <iostream> typedef double V;double finish(V x); double fnewton(V x); void divide(V a1,V a2); void newton(V a1); void cord(V a1,V a2); void newton2(V a1); using namespace std;double f(V x) retu

9、rn pow(x,5)-6*pow(x,3)+x-1; double fnewton(V x) if(5*pow(x,4)-9*pow(x,2)+1)!=0) return (pow(x,5)-3*pow(x,3)-1)/(5*pow(x,4)-9*pow(x,2)+1); elsereturn 1;void divide(V a1,V a2) for(;fabs(fabs(a1)-fabs(a2)>=1e-6;) if(f(a1)*f(a1+a2)/2)<0) a2=(a1+a2)/2;else a1=(a1+a2)/2; cout<<endl<<&quo

10、t;经过二分法求得的结果是："<<(a1+a2)/2<<endl<<endl; void newton(V a1) for(;fabs(fabs(a1)-fabs(a1-fnewton(a1)>=1e-6;) a1-=fnewton(a1); cout<<"经过Newton法得的结果是："<<a1<<endl<<endl;void cordline(V a1,V a2) int n=0; for(;fabs(f(a2-f(a2)*(a2-a1)/(f(a2)-f(a1)-0)

11、>=1e-6;n+) if(n=50) break; else if(f(a2-f(a2)*(a2-a1)/(f(a2)-f(a1)*f(a1)>0|(f(a2)-f(a1)!=0) a1=a2-f(a2)*(a2-a1)/(f(a2)-f(a1); else a2=a2-f(a2)*(a2-a1)/(f(a2)-f(a1); if(n=50) cout<<"使用双点弦法求解失败！"<<endl<<endl; elsecout<<"经过双点弦法求得的结果是："<<( a2-f(a2)

12、*(a2-a1)/(f(a2)-f(a1) )<<endl<<endl;void newton2(V a1) double t=1.0; for(;fabs(f(a1)<=fabs(f(a1-t*fnewton(a1);t/=2.0) for(;fabs(fabs(a1)-fabs(a1-t*fnewton(a1)>=1e-6;) a1-=t*fnewton(a1); cout<<"使用Newton下山法求得的结果是："<<a1<<endl<<endl;int main() double x

13、1,x2; cout<<"方程是：x5-3*x3+x-1= 0n请给出求根区间的左值和右值：n左值：" cin>>x1; cout<<"右值：" cin>>x2; if(f(x1)*f(x2)>0) cout<<"此区间上方程无根"<<endl; system("PAUSE"); return 0; divide(x1,x2); newton(x1+x2)/2); cordline(x1,x2); newton2(x1+x2)/2); s

14、ystem("PAUSE"); return 0; 实验2 线性方程组数值解-直接法一、实验目的：掌握Gauss消元法、列主元的Gauss消元法、Gauss-Jordan消元法、LU分解、改良平方根法等求解线性方程组的数值解的直接求解。二、实验内容：1. 用Gauss消元法、选主元的Gauss消元法求解以下线性方程组的解Ax=b，要求输出消元过程中系数矩阵的变化情况2. 用Gauss-Jordan消元法解题1中的线性方程组，并求A的逆矩阵A-13. 对上题中的线性方程组的系数矩阵A做LUDoolittle分解，要求输出L、U矩阵，并求线性方程组的解 (1) 高斯消元法算法分

15、析摘抄自课本在线性代数里就有用Gauss消元法线性方程组的解。通常，对线性方程组的增广矩阵A进行初等行变换化为上三角矩阵的方法，就成为Gauss消元法。对于任意一矩阵，均可以用Gauss消元法求解对应线性方程组的解。当用k表示消元过程的次序时，Gauss消元法的计算步骤为：1Gauss消元过程：设对计算2回代求解过程：(2) 列主元的Gauss消元法 算法分析列主元的Gauss消元法是在Gauss消元法的根底之上加上了列选主元的步2出绝对值最大的，然后通过行交换将其交换到的位置上。设主元在第个方程，即。假设，将和方程互易位置，使新的成为主元，然后继续进行，这一步骤称为列选主元。2.源代码#in

16、clude<iostream>#include<iomanip>using namespace std;void DisplayMatrix(double *a,int size)cout<<endl<<"="<<endl;for(int i=0;i<size;i+)for(int j=0;j<size+1;j+)cout<<setw(15)<<setprecision(8)<<aij;cout<<endl;cout<<endl<<

17、"="<<endl;void Gauss(double *a,int size)int i,j,k;double tmp;cout<<"方程组对应的增广矩阵为:"<<endl;DisplayMatrix(a,size);/输出矩阵,需自己实现该函数 /*消元,第k列对角线以下全消成0*/for(k=0;k<size-1;k+)for(i=k+1;i<size;i+)tmp=aik/akk;for(j=k;j<=size;j+)aij-=tmp*akj;cout<<endl<<&

18、quot;第"<<k+1<<"列消元后的矩阵为："<<endl;DisplayMatrix(a,size);/*回代求解*/for(i=size-1;i>=0;i-)for(j=size-1;j>i;j-)aisize-=ajsize*aij;aisize/=aii;/*输出方程组的解，放在数组最后一列中*/cout<<endl<<"方程的根为："<<endl;for(i=0;i<size;i+)cout<<"x"<&

19、lt;i+1<<"="<<setiosflags(ios:left)<<setw(15)<<setprecision(8)<<aisize;cout<<endl;int main() int n,i,j; cout<<"input the size of equations!"<<endl; cin>>n; double *p=new double *n; for( i=0;i<n;i+) pi=new double n+1; cout<

20、;<"Input the matrix of Equations!"<<endl; for(i=0;i<n;i+) for(j=0;j<=n;j+) cin>>pij; Gauss(p,n); for ( i=0;i<n;i+) delete pi; delete p; system("pause");截图：实验3 线性方程组数值解-迭代法一、实验目的：掌握Jaccobi迭代法、Guass-Sidel迭代法、松弛法求解线性方程组的数值解。二、实验内容：1分别用雅格比法与高斯赛德尔迭代法解以下方程组Ax=b，

21、研究其收敛性，上机验证理论分析是否正确，比拟它们的收敛速度，观察右端项对迭代收敛有无影响。(1)A行分别为A1=6,2,-1,A2=1,4,-2,A3=-3,1,4； b1=-3,2,4T, b2=100,-200,345T,(2) A行分别为A1=1,0,8,0.8,A2 0.8,1,0.8,A3=0.8,0.8,1；b1=3,2,1 T, b2=5,0,-10T,(3)A行分别为A1=1,3,A2=-7,1；b=4,6T,2. 选作松弛因子对SOR法收敛速度的影响。用SOR法求解方程组Ax=b,其中要求程序中不存系数矩阵A,分别对不同的阶数取进行迭代，记录近似解x(k)到达|x(k)-x(

22、k-1)|<10-6时所用的迭代次数k，观察松弛因子对收敛速度的影响，并观察当w£0或w³2会有什么影响？1Jacobi迭代法摘抄自课本 算法分析 线性方程组Ax=b，系数矩阵A非奇异，且aii0。即 可以简写为从上式中别离出变量，将它写成据此便建立了Jacobi迭代公式即可写为这就是解线性方程组的Jacobi迭代公式。也是Jacobi迭代法求解线性方程组的理论根底。源代码#include <iostream>#include<cmath>#include<iomanip>using namespace std;void jacco

23、bi(double *a,int size)cout<<"雅可比迭代法"double x100=0;double x_tp100=0;for(int i=0;i<17;i+)for(int i=0;i<size;i+) double tmp=aisize;for(int j=0;j<i;j+)tmp-=aij*xj;for(int j=i+1;j<size;j+)tmp-=aij*xj;x_tpi=tmp/aii;for(int i=0;i<size;i+) xi=x_tpi; cout<<endl<<&qu

24、ot;第"<<i+1<<"迭代：" for(int i=0;i<size;i+) cout<<"x"<<i+1<<"="<<left<<setw(9)<<setfill(' ')<<xi<<"|" cout<<endl; int main()int n,i,j;cout<<"请输入维数："cin>>n;doub

25、le *mat=new double *n;for( i=0;i<n;i+) mati=new double n+1; cout<<"请输入增广矩阵："<<endl; for(i=0;i<n;i+) for(j=0;j<=n;j+) cin>>matij; jaccobi(mat,n);gs(mat,n);for ( i=0;i<n;i+) delete mati; delete mat; system("PAUSE");return 0;2Guass_Sidel迭代法(摘抄自课本)算法分析在G

26、uass_Sidel迭代时，每次迭代充分利用当前最新的迭代值。在迭代收敛时，因新值比老值更准确些，求出新值后，用代替前一次的迭代值继续进行计算，这就充分利用新值建立起Guass_Sidel迭代公式。对于一般形式的方程组，Guass_Sidel迭代公式为源代码#include <iostream>#include<cmath>#include<iomanip>using namespace std;void gs(double *a,int size)cout<<"高斯赛德尔迭代法"double x100=0;double x_

27、tp100=0;for(int k=0;k<17;k+)for(int i=0;i<size;i+)double tmp=aisize;for(int j=0;j<i;j+)tmp-=aij*xj;for(int j=i+1;j<size;j+)tmp-=aij*xj;xi=tmp/aii;/for(int i=0;i<size;i+)/xi=x_tpi;cout<<"第"<<right<<setw(2)<<i+1<<"迭代："for(int i=0;i<si

28、ze;i+)cout<<"x"<<i+1<<"="<<left<<setw(9)<<setfill(' ')<<xi<<"|"cout<<endl;int main()int n,i,j;cout<<"请输入维数："cin>>n;double *mat=new double *n;for( i=0;i<n;i+)mati=new double n+1; cout&

29、lt;<"请输入增广矩阵："<<endl;for(i=0;i<n;i+)for(j=0;j<=n;j+)cin>>matij;jaccobi(mat,n);gs(mat,n);for ( i=0;i<n;i+)delete mati;delete mat;system("PAUSE");return 0;雅各布解的迭代过程a:高斯赛德尔迭代过程a：用雅各比经过 19 步迭代得到解为，。用高斯-塞德尔算法经过 11 步迭代得到的解为雅各布解的迭代过程b高斯赛德尔迭代过程b用雅各布经过 24 步迭代得到方程组的

30、解是x1 =36.37, x2 =-2.071, x3 =114.谱半径 (B) = 0.54用高斯-塞德尔经过 16 迭代得到方程组的解是x12高斯赛德尔迭代过程a,b用雅各布算法求的谱半径(B) = 1.6 > 1,由迭代矩阵的收敛条件知其不收敛。 用高斯-塞德尔经过 34 步迭代得到方程组的解是x1 =5.7691, x2 =0.7692, x3 =-4.2307 谱半径 (B) = 0.72 知用高斯-塞德尔法是收敛的。用雅各布用雅各布算法求的谱半径(B) = 1.6 > 1,由迭代矩阵的收敛条件知其不收敛。7用高斯-塞德尔经过 40 步迭代解方程组的解为x1谱半径 (B)

31、 = 0.72 得用高斯-塞德尔法是收敛的。(3)用雅各比算法迭代得到(B) = 4.6 >1,故迭代过程不收敛。 用高斯-塞德尔算法迭代得(B) = 21>1,迭代过程不收敛。 (1)在同样的精度要求和初始向量情况下,高斯-塞德尔的收敛速度要快于雅各比算法。 (2)迭代矩阵的收敛与否只和迭代矩阵B 的谱半径有关,同初始向量和方程组的右端项无关 。(3)一般情况下高斯-塞德尔收敛比雅各比快,但是也存在雅各比收敛而高斯-塞德尔不收 敛以及两种算法均不收敛的情况。实验4 Lagrange插值与三次样条插值效果的比拟一、实验目的：掌握Lagrange插值、三次样条插值的原理，了解Rung

32、e现象及防止方法二、实验内容：实际问题中两次过程记录测得(x,y)数据点如下：其中yi=f(xi) 第一次的数据为：x=-5.0000 -4.0000 -3.0000 ;y=-0.0000 2 -0.0004 -0.0366 -0.3679 0 0.3679 2第二次的数据为： x= ;y =1203 321针对第一次采集的数据请分别采用Lagrange插值与三次样条插值，估算在-5,5区间上各20等分点处的函数值，并将计算结果与第二次采集的结果进行比对，并分析误差原因。要求：以表格形式输出各计算结果，如：xyLarange(x)S(x)三、算法与步骤：1、拉格朗日插值：1输入数据点总数为n+

33、1即输入n的值，节点Xi和对应的Yi的值i从0到n，令Lnx= 0；2当i从0到n时，计算Li(x)= (X-Xj)/(Xi-Xj)其中j从0到n且j不等于i；令s=1，j从0到n：如果j=i，那么s=s；否那么令s=s*(x-xj/(xi-xj);循环结束令Lix=s，Ln(x) = (Ln(x)+Li(x)*Yi； 3)将待估侧值代入Lnx求值；2、三次样条插值假定有n+1个数据节点a. 计算步长 (i = 0, 1, , n-1)b. 将数据节点和指定的首位端点条件带入矩阵方程c. 解矩阵方程，求得二次微分值。该矩阵为三对角矩阵，具体求法参见我的上篇文章：三对角矩阵的求解。d.

34、 计算样条曲线的系数：其中i = 0, 1, , n-1e. 在每个子区间 中，创立方程 四、实验截图：1、拉格朗日插值2、三次样条插值附录：源程序清单1、拉格朗日插值源码#include<iostream>#include<iomanip>#define N 21 /插值节点数目using namespace std;void main()float aN; /所求点横坐标，将插值节点的横坐标复制float fx = 0, tmp = 1;int i, j;/差值节点横坐标float xN = -5.0000, -4.5000, -4.0000,

35、-3.5000, -3.0000,-2.5000, -2.0000, -1.5000, -1.0000, -0.5000,0, 0.5000, 1.0000, 1.5000, 2.0000, 4.0000, 4.5000,5.0000;/差值节点纵坐标float yN = -0.0000, -0.0001, -0.0002, -0.0003, -0.0004, -0.0048, -0.0366, , 0, 0.3894, 0.3679, 0.1581, , 0.0004, 0.0003, 0.0002, 0.0001, 0.0000 ;cout << "插值节点横坐标:n

36、"for (int i = 0; i < N; i+)cout << xi << " "cout << "n"cout << "插值节点纵坐标：n"for (int i = 0; i < N; i+)cout << yi << " "cout << "n"cout << "拉格朗日插值的结果为：" << endl;cout << &qu

37、ot;x"<< setw(20) << " y" << setw(20) << " Larange(x)" << endl;for (int i = 0; i < N; i+)ai = xi;for (int k = 0; k < N; k+)for (i = 0; i < N; i+)tmp = 1; /开始前复原tmpfor (j = 0; j < N; j+)if (i != j)tmp = tmp*(ak-xj) / (xi - xj);fx = fx

38、+ tmp*yi;cout <<xk<<setw(20)<<yk<<setw(20)<< fx << endl;fx = 0;2、三次样条插值源码#include<iostream>#include<iomanip>#include<cmath>#define N 21 /插值节点数目using namespace std; int main()/差值节点横坐标double xN = -5.0000, -4.5000, -4.0000, -3.5000, -3.0000,-2.5000,

39、 -2.0000, -1.5000, -1.0000, -0.5000,0, 0.5000, 1.0000, 1.5000, 2.0000, 4.0000, 4.5000,5.0000;/差值节点纵坐标double yN = -0.0000, -0.0001, -0.0002, -0.0003, -0.0004, -0.0048, -0.0366, , 0, 0.3894, 0.3679, 0.1581, , 0.0004, 0.0003, 0.0002, 0.0001, 0.0000 ;int choice = 0;double *xx,*a, *b, *a1, *b1, *h, *m;xx = new doubleN;a = new doubleN;b = new doubleN;a1 = new doubleN;b1 = new doubleN;h = new doubleN-1;m = new doubleN+1;for (int j = 0; j < N - 1; j+)hj = xj + 1 - xj;cout << "三次样条插值结果为：n