计算方法上机报告_董晓壮
(完整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 。
计算方法上机实验
1.拉格朗日插值多项式，用于离散数据的拟合#include <stdio.h>#include <conio.h>#include <alloc.h>float lagrange(float *x,float *y,float xx,int n) /*拉格朗日插值算法*/{ int i,j;float *a,yy=0.0; /*a作为临时变量，记录拉格朗日插值多项式*/a=(float *)malloc(n*sizeof(float));for(i=0;i<=n-1;i++){ a[i]=y[i];for(j=0;j<=n-1;j++)if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]);yy+=a[i];}free(a);return yy;}main(){ int i,n;float x[20],y[20],xx,yy;printf("Input n:");scanf("%d",&n);if(n>=20) {printf("Error!The value of n must in (0,20)."); getch();return 1;} if(n<=0) {printf("Error! The value of n must in (0,20)."); getch(); return 1;} for(i=0;i<=n-1;i++){ printf("x[%d]:",i);scanf("%f",&x[i]);}printf("\n");for(i=0;i<=n-1;i++){ printf("y[%d]:",i);scanf("%f",&y[i]);}printf("\n");printf("Input xx:");scanf("%f",&xx);yy=lagrange(x,y,xx,n);printf("x=%f,y=%f\n",xx,yy);getch();}2.牛顿插值多项式，用于离散数据的拟合#include <stdio.h>#include <conio.h>#include <alloc.h>void difference(float *x,float *y,int n){ float *f;int k,i;f=(float *)malloc(n*sizeof(float));for(k=1;k<=n;k++){ f[0]=y[k];for(i=0;i<k;i++)f[i+1]=(f[i]-y[i])/(x[k]-x[i]);y[k]=f[k];}return;}main(){ int i,n;float x[20],y[20],xx,yy;printf("Input n:");scanf("%d",&n);if(n>=20) {printf("Error! The value of n must in (0,20)."); getch(); return 1;} if(n<=0) {printf("Error! The value of n must in (0,20).");getch(); return 1;} for(i=0;i<=n-1;i++){ printf("x[%d]:",i);scanf("%f",&x[i]);}printf("\n");for(i=0;i<=n-1;i++){ printf("y[%d]:",i);scanf("%f",&y[i]);}printf("\n");difference(x,(float *)y,n);printf("Input xx:");scanf("%f",&xx);yy=y[20];for(i=n-1;i>=0;i--) yy=yy*(xx-x[i])+y[i];printf("NewtonInter(%f)=%f",xx,yy);getch();}3.高斯列主元消去法,求解其次线性方程组第一种#include<stdio.h>#include <math.h>#define N 20int main(){ int n,i,j,k;int mi,tmp,mx;float a[N][N],b[N],x[N];printf("\nInput n:");scanf("%d",&n);if(n>N){ printf("The input n should in(0,N)!\n");getch();return 1;}if(n<=0){ printf("The input n should in(0,N)!\n");getch();return 1;}printf("Now input a(i,j),i,j=0...%d:\n",n-1); for(i=0;i<n;i++){ for(j=0;j<n;j++)scanf("%f",&a[i][j]);}printf("Now input b(i),i,j=0...%d:\n",n-1);for(i=0;i<n;i++)scanf("%f",&b[i]);for(i=0;i<n-2;i++){ for(j=i+1,mi=i,mx=fabs(a[i][j]);j<n-1;j++) if(fabs(a[j][i])>mx){ mi=j;mx=fabs(a[j][i]);}if(i<mi){ tmp=b[i];b[i]=b[mi];b[mi]=tmp;for(j=i;j<n;j++){ tmp=a[i][j];a[i][j]=a[mi][j];a[mi][j]=tmp;}}for(j=i+1;j<n;j++){ tmp=-a[j][i]/a[i][i];b[j]+=b[i]*tmp;for(k=i;k<n;k++)a[j][k]+=a[i][k]*tmp;}}x[n-1]=b[n-1]/a[n-1][n-1];for(i=n-2;i>=0;i--){ x[i]=b[i];for(j=i+1;j<n;j++)x[i]-=a[i][j]*x[j];x[i]/=a[i][i];}for(i=0;i<n;i++)printf("Answer:\n x[%d]=%f\n",i,x[i]); getch();return 0;}第二种#include<math.h>#include<stdio.h>#define NUMBER 20#define Esc 0x1b#define Enter 0x0dfloat A[NUMBER][NUMBER+1] ,ark;int flag,n;exchange(int r,int k);float max(int k);message();main(){float x[NUMBER];int r,k,i,j;char celect;clrscr();printf("\n\nUse Gauss.");printf("\n\n1.Jie please press Enter."); printf("\n\n2.Exit press Esc.");celect=getch();if(celect==Esc)exit(0);printf("\n\n input n=");scanf("%d",&n);printf(" \n\nInput matrix A and B:"); for(i=1;i<=n;i++){printf("\n\nInput a%d1--a%d%d and b%d:",i,i,n,i);for(j=1;j<=n+1;j++) scanf("%f",&A[i][j]);}for(k=1;k<=n-1;k++){ark=max(k);if(ark==0){printf("\n\nIt's wrong!");message();}else if(flag!=k)exchange(flag,k);for(i=k+1;i<=n;i++)for(j=k+1;j<=n+1;j++)A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k];}x[n]=A[n][n+1]/A[n][n];for( k=n-1;k>=1;k--){float me=0;for(j=k+1;j<=n;j++){me=me+A[k][j]*x[j];}x[k]=(A[k][n+1]-me)/A[k][k];}for(i=1;i<=n;i++){printf(" \n\nx%d=%f",i,x[i]);}message();}exchange(int r,int k){int i;for(i=1;i<=n+1;i++)A[0][i]=A[r][i];for(i=1;i<=n+1;i++)A[r][i]=A[k][i];for(i=1;i<=n+1;i++)A[k][i]=A[0][i];}float max(int k){int i;float temp=0;for(i=k;i<=n;i++)if(fabs(A[i][k])>temp){temp=fabs(A[i][k]);flag=i;}return temp;}message(){printf("\n\n Go on Enter ,Exit press Esc!");switch(getch()){case Enter: main();case Esc: exit(0);default:{printf("\n\nInput error!");message();} }}4.龙贝格求积公式，求解定积分#include<stdio.h>#include<math.h>#define f(x) (sin(x)/x)#define N 20#define MAX 20#define a 2#define b 4#define e 0.00001float LBG(float p,float q,int n){ int i;float sum=0,h=(q-p)/n;for (i=1;i<n;i++)sum+=f(p+i*h);sum+=(f(p)+f(q))/2;return(h*sum);}void main(){ int i;int n=N,m=0;float T[MAX+1][2];T[0][1]=LBG(a,b,n);n*=2;for(m=1;m<MAX;m++){ for(i=0;i<m;i++)T[i][0]=T[i][1];T[0][1]=LBG(a,b,n);n*=2;for(i=1;i<=m;i++)T[i][1]=T[i-1][1]+(T[i-1][1]-T[i-1][0])/(pow(2,2*m)-1);if((T[m-1][1]<T[m][1]+e)&&(T[m-1][1]>T[m][1]-e)){ printf("Answer=%f\n",T[m][1]); getch();return ;}}}5.牛顿迭代公式，求方程的近似解#include<stdio.h>#include<math.h>#include<conio.h>#define N 100#define PS 1e-5#define TA 1e-5float Newton(float (*f)(float),float(*f1)(float),float x0 ) { float x1,d=0;int k=0;do{ x1= x0-f(x0)/f1(x0);if((k++>N)||(fabs(f1(x1))<PS)){ printf("\nFailed!");getch();exit();}d=(fabs(x1)<1?x1-x0:(x1-x0)/x1);x0=x1;printf("x(%d)=%f\n",k,x0);}while((fabs(d))>PS&&fabs(f(x1))>TA) ;return x1;}float f(float x){ return x*x*x+x*x-3*x-3; }float f1(float x){ return 3.0*x*x+2*x-3; }void main(){ float f(float);float f1(float);float x0,y0;printf("Input x0: ");scanf("%f",&x0);printf("x(0)=%f\n",x0);y0=Newton(f,f1,x0);printf("\nThe root is x=%f\n",y0); getch();}。
数值计算方法上机实验报告
数值计算方法上机实验报告
一、实验目的
本次实验的主要目的是熟悉和掌握数值计算方法，学习梯度下降法的
原理和实际应用，熟悉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库的使用方法。
计算方法上上机实习报告
计算方法上上机实习报告在本次计算方法的上机实习中，我深入体验了数值计算的魅力和挑战，通过实际操作和实践，对计算方法有了更深刻的理解和认识。
实习的目的在于将课堂上学到的理论知识运用到实际的计算中，熟悉各种数值算法的实现过程，提高编程能力和解决实际问题的能力。
我们使用了具体编程语言和软件名称进行编程和计算。
在实习过程中，我首先接触到的是数值逼近的相关内容。
通过多项式插值和曲线拟合的练习，我明白了如何用简单的函数去近似复杂的曲线。
例如，拉格朗日插值法和牛顿插值法让我能够根据给定的离散数据点构建出一个连续的函数，从而对未知点进行预测。
在实际操作中，我需要仔细处理数据的输入和输出，以及算法中的细节，如边界条件和误差控制。
数值积分是另一个重要的部分。
通过梯形公式和辛普森公式，我学会了如何对给定的函数进行数值积分。
在编程实现时，要合理地选择积分区间和步长，以达到所需的精度。
同时，我也了解到了数值积分方法的误差来源和误差估计方法，这对于评估计算结果的可靠性非常重要。
线性方程组的求解是计算方法中的核心内容之一。
我分别使用了高斯消元法和迭代法（如雅克比迭代法和高斯赛德尔迭代法）来求解线性方程组。
在实际编程中，我深刻体会到了算法的效率和稳定性的重要性。
对于大规模的线性方程组，选择合适的算法可以大大提高计算速度和精度。
在非线性方程求根方面，我运用了二分法、牛顿法和割线法等方法。
这些方法各有特点，二分法简单但收敛速度较慢，牛顿法收敛速度快但需要计算导数。
在实际应用中，需要根据方程的特点和求解的要求选择合适的方法。
在实习中，我也遇到了不少问题和挑战。
首先是编程中的错误，如语法错误、逻辑错误等，这需要我耐心地调试和修改代码。
其次，对于一些复杂的算法，理解其原理和实现细节并不容易，需要反复查阅资料和思考。
还有就是数值计算中的误差问题，有时候由于误差的积累，导致计算结果与预期相差较大，需要通过调整算法参数或者采用更精确的算法来解决。
《数值计算方法》上机实验报告
《数值计算方法》上机实验报告华北电力大学实验名称数值il•算方法》上机实验课程名称数值计算方法专业班级：电力实08学生姓名：李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期：2010年04月华北电力大学实验报告数值计算方法上机实验报吿一.各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程*对于非线性方程，若已知根的一个近似值，将在处展开成一阶xxfx ()0, fx ()xkk泰勒公式"f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2!忽略高次项，有，fxfxfxxx 0 ()()()，，, kkk右端是直线方程，用这个直线方程来近似非线性方程。
将非线性方程的**根代入，即fx ()0, X ，* fxfxxx 0 0 0 0, ,, kkkfx 0 fx 0 0,解出fX 0 *k XX,, k' fx 0 k水将右端取为，则是比更接近于的近似值，即xxxxk, Ik, Ikfx ()k 八XX, Ikk* fx()k这就是牛顿迭代公式。
,2,计算机程序框图:，见,,3,输入变量、输出变量说明:X输入变量:迭代初值，迭代精度，迭代最大次数，\0输出变量:当前迭代次数，当前迭代值xkl,4,具体算例及求解结果:2/16华北电力大学实验报吿开始读入l>k/fx()0?,0fx 0 Oxx，，01* fx ()0XX,,，?10kk, ,1，kN, ?xx, 10输出迭代输出X输出奇异标志1失败标志,3,输入变量、输出变量说明: 结束例：导出计算的牛顿迭代公式，并il •算。
（课本P39例2-16） 115cc （0）, 求解结果:10. 75000010.72383710. 72380510. 7238052、列主元素消去法求解线性方程组，1,算法原理:高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘 -个 方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上 对上三角3/16华北电力大学实验报告方程组求解。
计算方法与实习上机实验报告
计算方法与实习上机实验报告一、引言本文旨在介绍和展示我们在“计算方法”课程中的实习上机实验环节所完成的一些关键任务和所取得的成果。
该实验课程的目标是让我们更深入地理解和应用各种计算方法，并在实际操作中提高我们的编程和问题解决能力。
二、实验内容与目标实验的主要内容是利用各种计算方法解决实际数学问题。
我们被要求使用编程语言（如Python或Java）来实现和解决这些问题。
这些问题包括使用牛顿法求解平方根，使用蒙特卡洛方法计算圆周率，以及使用最优化方法求解函数的最小值等。
实验的目标不仅是让我们掌握计算方法的基本理论，更是要让我们能够在实际操作中运用这些方法。
我们需要在实习过程中，通过与同伴们合作，共同解决问题，提高我们的团队合作能力和问题解决能力。
三、实验过程与问题解决策略在实验过程中，我们遇到了许多问题，如编程错误、理解困难和时间压力等。
我们通过相互讨论、查阅资料和寻求教师帮助等方式，成功地解决了这些问题。
例如，在实现牛顿法求解平方根时，我们一开始对导数的计算和理解出现了一些错误。
但我们通过查阅相关资料和讨论，最终理解了导数的正确计算方法，并成功地实现了牛顿法。
四、实验结果与结论通过这次实习上机实验，我们不仅深入理解了计算方法的基本理论，还在实际操作中提高了我们的编程和问题解决能力。
我们的成果包括编写出了能有效求解平方根、计算圆周率和求解函数最小值的程序。
这次实习上机实验非常成功。
我们的团队不仅在理论学习和实践操作上取得了显著的进步，还在团队合作和问题解决方面积累了宝贵的经验。
这次实验使我们对计算方法有了更深的理解和认识，也提高了我们的编程技能和解决问题的能力。
五、反思与展望回顾这次实验，我们意识到在实验过程中，我们需要更好地管理我们的时间和压力。
在解决问题时，我们需要更有效地利用我们的知识和资源。
在未来，我们希望能够更加熟练地运用计算方法，并能够更有效地解决问题。
我们也希望能够将所学的计算方法应用到更广泛的领域中，如数据分析、科学研究和工业生产等。
计算方法与实习上机报告
计算方法与实习——上机报告学院：电子工程学院学号：姓名：刘波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.逆序计算结果分析：两种不同顺序计算结果是一样的，顺序计算误差从一开始就很小，而逆序计算误差最开始十分大，后来结果正确。
计算机上机报告
计算方法上机实验报告上课时间: 2014-2015学年秋学期, 6～14周一. 拉格朗日插值------------------------------------------------------1二. 牛顿插值------------------------------------------------------------3三. 改进欧拉法---------------------------------------------------------5四. 四阶龙格-库塔-----------------------------------------------------7五. 牛顿迭代------------------------------------------------------------9 六．复化Simpson公式------------------------------------------------11七. Romberg算法------------------------------------------------------14八. Seidel迭代法------------------------------------------------------17九. Gauss列主元消去法----------------------------------------------20一.拉格朗日插值1.程序代码#include<iostream.h>void Lagrange(){int i=0;double a[10],b[10],L,L1,L2,L3,L4,x;cout<<"x=";for(i=0;i<4;i++){cin>>a[i];}cout<<"y=";for(i=0;i<4;i++){cin>>b[i];}cout<<"x=";cin>>x;L1=(x-a[1])*(x-a[2])*(x-a[3])*b[0]/(a[0]-a[1])/(a[0]-a[2])/(a[0]-a[3]);L2=(x-a[0])*(x-a[2])*(x-a[3])*b[1]/(a[1]-a[0])/(a[1]-a[2])/(a[1]-a[3]);L3=(x-a[0])*(x-a[1])*(x-a[3])*b[2]/(a[2]-a[0])/(a[2]-a[1])/(a[2]-a[3]);L4=(x-a[0])*(x-a[1])*(x-a[2])*b[3]/(a[3]-a[0])/(a[3]-a[1])/(a[3]-a[2]);L=L1+L2+L3+L4;cout<<"L="<<L;}void main(){Lagrange();cout<<endl;}2.例子3.运行结果二. 牛顿插值1.程序代码#include <iostream.h>#include<conio.h>void main(){int n,i,j;double A[50][50],*x,*y;cout<<"请输入插值节点数: ";cin>>n;x=new double[n];y=new double[n];cout<<"请输入这"<<n<<"个插值节点（xi,yi）: "<<endl;for(i=0;i<=n-1;i++)cin>>x[i]>>y[i];double K=1,xx,N=0,P;for(i=0;i<=n-1;i++){A[i][0]=x[i];A[i][1]=y[i];}for(j=2;j<=n;j++)for(i=1;i<=n-1;i++){A[i][j]=(A[i][j-1]-A[i-1][j-1])/(A[i][0]-A[i-j+1][0]);}for(i=0;i<=n-1;i++)cout<<"输出第"<<i+1<<"阶差商为: "<<A[i][i+1]<<endl;cout<<"请输入预求值x=";cin>>xx;for(i=0;i<n-1;i++){K*=xx-x[i];N+=A[i+1][i+2]*K;P=A[0][1]+N;}cout<<"Newton插值结果为: y="<<P<<endl;getch();}2.例子3.运行结果三. 改进欧拉法1.程序代码#include<iostream.h>#include<conio.h>double fun(double x,double y){return(-0.9*y/(1+2*x));}void main(){double a,b,*y,h,*x,yp,yc;int n,k;cout<<"常微分方程为y'=-0.9*y/(1+2*x)"<<endl;cout<<"其中0<=x<=1"<<endl;cout<<"初值为y(0)=1"<<endl;cout<<"请输入计算区间(a,b): ";cin>>a>>b;cout<<"请输入步长h: ";cin>>h;cout<<"请输入计算次数: ";cin>>n;y=new double[n];x=new double[n];cout<<"请输入初值y(0)=";cin>>y[0];x[0]=a;for(k=0;k<=n;k++){yp=y[k]+h*fun(x[k],y[k]);yc=y[k]+h*fun(x[k]+h,yp);y[k+1]=0.5*(yp+yc);x[k+1]=x[k]+h;}cout<<"迭代结果为: "<<endl;for(k=0;k<=n;k++)cout<<"x"<<k<<"="<<x[k]<<'\t'<<"y"<<k<<"="<<y[k]<<endl;getch();}2.例子3.运行结果四. 四阶龙格-库塔1.程序代码#include<iostream.h>#include<conio.h>double fun(double x,double y){return(x-y);}void main(){double a,b,*y,h,x,k1,k2,k3,k4;int n,k;cout<<"常微分方程为y'=x-y"<<endl;cout<<"其中0<=x<=1"<<endl;cout<<"初值为y(0)=0"<<endl;cout<<"请输入计算区间(a,b): ";cin>>a>>b;cout<<"请输入步长h: ";cin>>h;cout<<"请输入计算次数: ";cin>>n;y=new double[n];cout<<"请输入初值y(0): ";cin>>y[0];x=a;cout<<"Runge-Kutta迭代法结果为: "<<endl;cout<<"x0="<<x<<'\t'<<"y0="<<y[0]<<endl;for(k=0;k<=n-1;k++){k1=fun(x,y[k]);k2=fun(x+h/2,y[k]+k1*h/2);k3=fun(x+h/2,y[k]+k2*h/2);k4=fun(x+h,y[k]+k3*h);y[k+1]=y[k]+(h/6)*(k1+2*(k2+k3)+k4);cout<<"x"<<k+1<<"="<<x+h<<'\t'<<"y"<<k+1<<"="<<y[k+1]<<endl;x=x+h;}getch();}2.例子3.运行结果五. 牛顿迭代法1.程序代码(C++代码)#include<iostream>#include<cmath>using namespace std;double newtondiedai(double a,double b,double c,double d,double x); int main(){double a,b,c,d;double x=1.5;cout<<"请依次输入方程四个系数: ";cin>>a>>b>>c>>d;x=newtondiedai(a,b,c,d,x);cout<<x<<endl;return 0;}double newtondiedai(double a,double b,double c,double d,double x) {while(abs(a*x*x*x+b*x*x+c*x+d)>0.000001){x=x-(a*x*x*x+b*x*x+c*x+d)/(3*a*x*x+2*b*x+c);}return x;}2.例子3.运行结果六. 复化Simpson公式1.程序代码（C++代码）#include<iostream.h>#include<math.h>double function1(double x)//被积函数{double s;s=x/(4+x*x);return s;}double function2(double x)//被积函数{double s;s=sqrt(x);return s;}double ReiterationOfSimpson(double a,double b,double n,double f(double x))//复化Simpson公式{double h,fa,fb,xk,xj;h=(b-a)/n;fa=f(a);fb=f(b);double s1=0.0;double s2=0.0;for(int k=1;k<n;k++){xk=a+k*h;s1=s1+f(xk);}for(int j=0;j<n;j++){xj=a+(j+0.5)*h;s2=s2+f(xj);}double sn;//和sn=h/6*(fa+fb+2*s1+4*s2);//复化Simpson公式return sn;}main(){double a,b,Result,n;cout<<"请输入积分下限:"<<endl;cin>>a;cout<<"请输入积分上限:"<<endl;cin>>b;cout<<"请输入分割区间数n:"<<endl;cin>>n;cout<<"复化Simpson 公式计算结果:";Result=ReiterationOfSimpson(a, b, n,function1); cout<<Result<<endl;}2.例子 （620(n 3)4x I dx x ==+⎰）3.运行结果七. Romberg算法1.程序代码（C++代码）#include<iostream>#include<cmath>using namespace std;#define f(x) (4/(1+x*x))#define epsilon 0.0001#define MAXREPT 10double Romberg(double aa,double bb) { int m,n;double h,x;double s,q;double ep;double *y =new double[MAXREPT]; double p;h=bb-aa;y[0]=h*(f(aa)+f(bb))/2.0;m=1;n=1;ep=epsilon+1.0;while((ep>=epsilon)&&(m<MAXREPT)) { p =0.0;for(int i=0;i<n;i++){ x=aa+(i+0.5)*h;p=p+f(x);}p=(y[0] + h*p)/2.0;s=1.0;for(int k=1;k<=m;k++){ s=4.0*s;q=(s*p-y[k-1])/(s-1.0); y[k-1]=p;p=q;}p=fabs(q-y[m-1]);m=m+1;y[m-1]=q;n=n+n;h=h/2.0;}return (q);}int main(){double a,b;cout<<"Romberg积分,请输入积分范围a,b:"<<endl; cin>>a>>b;cout<<"积分结果:"<<Romberg(a,b)<<endl;system("pause");return 0;}2.例子（1241I dxx=+⎰）3.运行结果八. Seidel迭代法1.程序代码（C++代码）# include <math.h># include <stdio.h># define max 100# define EPS 1e-6float a[3][3]={{10,-1,-2},{-1,10,-2},{-1,-1,5}}; float b[3]={7.2,8.3,4.2};float x[3]={0,0,0};float y[3];float S(int m){int n;float S=0;float y;for(n=0;n<3;n++){if(m==n){}else{S+=a[m][n]*x[n];}}y=(b[m]-S)/a[m][m];return y;}void main(){int i;int F,T=1,k=1;do{F=0;if(T){T=0;}else{k++;for(i=0;i<3;i++){x[i]=y[i];}}for(i=0;i<3;i++){y[i]=S(i);}for(i=0;i<3;i++){if(fabs(x[i]-y[i])>EPS){F=1;}}printf("%d\n",k);}while(((F==1)&&k!=max));printf("迭代次数:%d\n",k);for(i=0;i<3;i++){printf("x[%d]=%1.6f\n",i,y[i]); }}2.例子（1012237.21102238.31253 4.2x x xx x xx x x--=⎧⎪-+-=⎨⎪--+=⎩）3.运行结果九. Gauss列主元消去法1.程序代码（C++代码）#include <stdio.h>#include<conio.h>#include <math.h>#define max_dimension 20int n;static float a[max_dimension][max_dimension]; static float b[max_dimension];static float x[max_dimension];void main(){int i;int j;int d;int row;float temp;float known_items;float l[max_dimension][max_dimension];printf("请输入阶数:");scanf("%d",&n);printf("\n");printf("请输入系数矩阵的值: ");printf("\n");for(i=0; i<n; i++){ printf("输入第%d行的值:",i+1);for (j=0; j<n; j++){scanf("%f",&a[i][j]);}printf("\n");}printf("请输入常数项的值: ");for(i=0; i<n; i++)scanf("%f",&b[i]);printf("\n");for(d=0; d<n-1; d++){row=d;for(i=d+1; i<n; i++){if(fabs(a[i][d])>fabs(a[row][d]))row=i;}if(row!=d){for(j=d; j<n; j++){temp=a[row][j];a[row][j]=a[d][j];a[d][j]=temp;}temp=b[row];b[row]=b[d];b[d]=temp;}for(i=d+1; i<n; i++){l[i][d]=-a[i][d]/a[d][d];for (j=d; j<n; j++){a[i][j]=a[i][j]+a[d][j]*l[i][d];}b[i]=b[i]+b[d]*l[i][d];}}printf("\n");for (i=n-1; i>-1; i--){known_items=0;for(j=1; j<n-i; j++){known_items=known_items+a[i][i+j]*x[i+j]; }x[i]=(b[i]-known_items)/a[i][i];}printf("方程组的根为:\n\n");for(i=0; i<n; i++)printf("%.5f ",x[i]);printf("\n\n");getch();}2.例子3.运行结果。
计算方法上机作业
计算方法上机报告姓名：学号：班级：上课班级:说明：本次上机实验使用的编程语言是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。
计算方法上机报告
计算方法上机报告硕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运行结果由图可知，三次样条插值拟合数据点比较好，没有龙格现象。
数值计算方法上机实验报告1
数值计算方法上机实验报告1华北电力大学上机实验报告课程名称：数值计算方法专业班级：学生姓名：学号：指导教师：张建成实验目的：复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。
利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。
上机练习任务：利用计算机基本C 语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。
掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。
一、列主元素消去法求解线性方程组 1、算法原理为避免绝对值很小的元素作为主元，在每次消元之前增加一个选主元的过程，将绝对值大的元素交换到主对角线的位置。
列主元素消元法是当变换到第k 步时，从k 列的kk a 及以下的各元素中选取绝对值最大的元素，然后通过二交换将其交换到kk a 的位置上。
2、输入输出变量ija ：为系数矩阵的各个系数K ：表示到第k 步消元 3、具体算例输入增广矩阵为： 3 1 2 -3 8 2 1 3 22 3 2 1 28解得：1x =6，2x =4，3x =2；二、LU 分解法求解线性方程组1、算法原理应用高斯消去法解n 阶线性方程Ax b =经过1n -步消去后得出一个等价的上三角形方程组()()n n A x b =，对上三角形方程组用逐步回代就可以求出解来。
这个过程也可通过矩阵分解来实现。
将非奇异阵分解成一个下三角阵L 和上三角阵U 的乘积A LU =称为对矩阵A 的三角分解，又称LU 分解。
根据LU 分解,将Ax b =分解为Ly bUx y =??=?形式,简化了求解问题。
2、输入输出变量ij a 为系数矩阵元素i b 为常数矩阵系数,i j i jl u 分别为下、上三角矩阵元素 k 表示第k 步消元 3、具体算例输入增广矩阵 3 2 3 4 39 3 -2 2 14 4 2 3 43 解得： 6 5 3三、拉格朗日插值1、算法原理设函数()y f x =在区间[a,b]上有节点01,,,,n x x x 上的函数值，构造一个次数不超过n次的代数多项式1110()n n n n p x a x a x a x a --=++++ ,使 (),0,1,,i i P x y i n == 。
计算方法实验上机报告(完整版)
简单迭代法#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];}。
计算方法上机实验报告
《计算方法》上机实验报告班级：XＸＸXXX小组成员:XXＸXXＸＸXXXＸXXXXXＸＸXXXＸXXXXXＸ任课教师：XXX二〇一八年五月二十五日前言通过进行多次得上机实验,我们结合课本上得内容以及老师对我们得指导,能够较为熟练地掌握Ｎewton迭代法、Ｊacｏｂi迭代法、Gausｓ－Seｉｄｅl 迭代法、Ｎewｔon 插值法、Ｌａgｒanｇe 插值法与Ｇaｕｓs 求积公式等六种算法得原理与使用方法，并参考课本例题进行了MATＬAB程序得编写。
以下为本次上机实验报告,按照实验内容共分为六部分.实验一：一、实验名称及题目:Newｔon迭代法例2、7(P3８）:应用Newｔon迭代法求在附近得数值解，并使其满足、二、解题思路：设就是得根,选取作为初始近似值，过点做曲线得切线,得方程为,求出与轴交点得横坐标,称为得一次近似值,过点做曲线得切线,求该切线与轴得横坐标称为得二次近似值，重复以上过程，得得近似值序列，把称为得次近似值，这种求解方法就就是牛顿迭代法。
三、Matlab程序代码：fｕncｔion newtｏn_itｅrａｔioｎ(x０,tol）syms z %定义自变量fｏrmaｔlｏnｇ％定义精度ｆ=ｚ*z＊z-z－1;f1=diff(f）；％求导y=suｂs(ｆ，ｚ,x0);y1=sｕbs(f１,ｚ，x0)；%向函数中代值x1=x0－y/y1; k=1；wｈilｅａｂｓ(ｘ1—x0)〉=ｔolx０=x1;y=subｓ（f,z,ｘ0)；y1=subｓ(ｆ1,ｚ,ｘ0);x1＝x0-ｙ/y1；k=k+1;ｅndｘ=douｂle（x１)Ｋ四、运行结果：实验二:一、实验名称及题目：Jacｏbｉ迭代法例３、7(P74）:试利用Jacobi迭代公式求解方程组要求数值解为方程组得精确解、二、解题思路：首先将方程组中得系数矩阵分解成三部分,即:，为对角阵,为下三角矩阵,为上三角矩阵。
之后确定迭代格式,,（, 即迭代次数),称为迭代矩阵。
计算方法(A)上机实验2015
计算方法(A) 上机实验报告一、共轭梯度法求解线型方程组 1.算法原理共轭梯度法求解关键在于搜索方向与步长，其思路是在点)(k X 处选取搜索方向)(k d , 使其与前一次的搜索方向)1(-k d 关于A 共轭，即(1)()(1),0k k k d d Ad --<>=然后从点)(k X 出发，沿方向)(k d 求得)(X f 的极小值点)1(+k X , 即)(min )()()(0)1(k d X f X f k k λλ+=>+)()1()1()()()()1(,,k k k k k k k d Ad d d AX b XX><>-<+=--+ 注意到)(k d 的选取不唯一，我们可取)1(1)()()(--+-∇=k k k k d X f d β由共轭的定义0,)1()(>=<-k k Add 可得：><><-=----)1()1()1()(1,,k k k k k Ad d Ad r β 共轭梯度法的计算过程如下： 第一步：取初始向量)0(X , 计算⎪⎪⎩⎪⎪⎨⎧+=><><-=-=-∇==(0)0(0)(1))0()0()0()0(0(0)(0)(0)(0)d X X ,,X )X (r d λλAd d Ad r A b f 第1+k 步：计算⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧+=><><-=+=><><-=-=-∇=+------(k)0(k)1)(k )()()()()1(1(k))()1()1()1()(1(k)(k)(k)d X X ,,r ,,X )X (r λλββk k k k k k k k k k k k k Ad d Ad r d d Ad d Adr A b f2. 程序框图3. 程序说明及算例程序为matlab 中的自定义函数文件，求解时调用即可，如C_gradient（[],x0,esp ）,其中矩阵为增广矩阵，迭代初值x0及精度esp 为选择性输入。
(完整word版)计算方法上机报告-备份
计算方法上机作业报告姓名：李小盼班级：计算方法B3班学号：6230489477专业：热能工程2016年《计算方法B》上机题目一.计算机语言要求使用语言不做限制，可以使用C、C++、FORTRAN、VC、VB、C#、Matlab、PHP、JavaScript等语言。
二.实习报告内容上机题目完成后，必须交一份上机报告。
上机报告中应对每道题目包括：（1）上机题目内容;（2）详细说明实现的思想、算法依据、算法实现的结构;（3）详细完整的源程序，并附相关的注释说明;（4）给出必要的计算结果（若数据量大，则只需列出主要的数据内容），并对结果进行分析;（5）对上机中出现的问题进行分析总结；三.实习报告要求1.提供一份完整的上机报告的电子文档；然后再提供一份与电子文档对应的纸质上机报告。
2.电子文档中提供上机过程中的所有源程序、输出数据，以及可以运行的文件。
如果程序的运行环境特殊，请附上运行程序所需要的软件环境。
3.上机报告严禁抄袭，如发现有抄袭现象，所有涉及抄袭的上机报告将被以作弊处理，并按零分处理,不再另行通知。
4.上机报告电子版一、二、三班发送到邮箱:xjtujsff@,上机作业纸面作业请送到:理科楼338。
上机实习题目1.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:米)如下表所示：(1)请用合适的曲线拟合所测数据点；(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；解：1.1 实现思想选用曲线拟合数据点时，一方面要满足插值条件，即保证所得曲线经过以上所有数据点，另一方面也要保证曲线具有足够的“光滑性”，故这里采用三次样条插值法构造插值函数。
为构成封闭方程组，边界条件选取自然三次样条，以获得三次样条插值的关键参数M0和M n。
1.2 算法的依据以及结构三次样条插值具有较好的光滑性，可以用于对数据点的光滑拟合。
哈工大计算方法上机实验报告
实验报告一题目： 非线性方程求解摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。
本实验采用两种常见的求解方法二分法和Newton 法及改进的Newton 法。
前言：（目的和意义）掌握二分法与Newton 法的基本原理和应用。
数学原理：对于一个非线性方程的数值解法很多。
在此介绍两种最常见的方法：二分法和Newton 法。
对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b ]上连续，f(a)f(b)<0，且f(x)在[a,b ]内仅有一个实根x *，取区间中点c ，若，则c 恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c ]和[c,b ]中的哪一个，从而得出新区间，仍称为[a,b ]。
重复运行计算，直至满足精度为止。
这就是二分法的计算思想。
Newton 法通常预先要给出一个猜测初值x 0，然后根据其迭代公式)()('1k k k k x f x f x x -=+ 产生逼近解x *的迭代数列{x k }，这就是Newton 法的思想。
当x 0接近x *时收敛很快，但是当x 0选择不好时，可能会发散，因此初值的选取很重要。
另外，若将该迭代公式改进为)()('1k k k k x f x f r x x -=+ 其中r 为要求的方程的根的重数，这就是改进的Newton 法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton 法快的多。
程序设计：本实验采用Matlab 的M 文件编写。
其中待求解的方程写成function 的方式，如下 function y=f(x);y=-x*x-sin(x);写成如上形式即可，下面给出主程序。
二分法源程序：clear%%%给定求解区间b=1.5;a=0;%%%误差R=1;k=0;%迭代次数初值while (R>5e-6) ;c=(a+b)/2;if f12(a)*f12(c)>0;a=c;elseb=c;endR=b-a;%求出误差k=k+1;endx=c%给出解Newton法及改进的Newton法源程序：clear%%%% 输入函数f=input('请输入需要求解函数>>','s')%%%求解f(x)的导数df=diff(f);%%%改进常数或重根数miu=2;%%%初始值x0x0=input('input initial value x0>>');k=0;%迭代次数max=100;%最大迭代次数R=eval(subs(f,'x0','x'));%求解f(x0)，以确定初值x0时否就是解while (abs(R)>1e-8)x1=x0-miu*eval(subs(f,'x0','x'))/eval(subs(df,'x0','x'));R=x1-x0;x0=x1;k=k+1;if (eval(subs(f,'x0','x'))<1e-10);breakendif k>max ;%如果迭代次数大于给定值，认为迭代不收敛，重新输入初值 ss=input('maybe result is error,choose a new x0,y/n?>>','s');if strcmp(ss,'y')x0=input('input initial value x0>>');k=0;elsebreakendendendk;%给出迭代次数x=x0；%给出解结果分析和讨论：1. 用二分法计算方程02sin 2=-x x 在[1，2]内的根。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算方法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） (),,P i j θ是正交矩阵，()()1T,,,,P i j P i j θθ−=；（2） 仅改变向量的第,i j 个元素； （3） 仅改变矩阵的第,i j 行元素；（4） 可以使用有限个吉文斯矩阵乘积P ，使得12Px x e =。
()111,,111n ncsP i j sc θ⨯⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪=⎪ ⎪ ⎪− ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭i j i j（1.1）矩阵A 按列分块表示为：()12,,,n a a a ，利用性质1,2,4，可以对每一列的各个元素进行计算消去，最终可以得到图1.1的QR 矩阵。
1.1.2 基于豪斯霍尔德变换的QR 分解给定一个n 阶的向量u ，满足21u=，矩阵T 2H I uu =−就是所说的豪斯霍尔德变换矩阵，也称为初等反射阵。
该矩阵存在一定的几何意义，当3n =时，如图1.2所示，对于空间给定的向量x 可以分解为12x x x =+，其中，u 为垂直于π平面的单位法向量，于是有：120x u x ku ==， （1.2）x 1x 2x uHx2x −图1.2 空间向量分解图将空间给定的向量x 和豪斯霍尔德矩阵相乘，有：()()()T 11T T22222Hx I uu x Hx I uu x Hx I uu x ⎧=−⎪=−⇒⎨=−⎪⎩ （1.3）结合式（1.2）可以进行化简，得到：1122Hx x Hx x =⎧⎨=−⎩ （1.4） 得到的向量关于平面对称，因此也称为初等反射矩阵。
豪斯霍尔德矩阵同样存在着优良的性质： （1） 矩阵H 为对称正交矩阵T 1H H H −==；（2） 对于任意给定的单位向量u ，存在初等反射矩阵H ，使得Hx u σ=； 矩阵A 按列分块表示为：()12,,,n a a a ，利用上述性质，可以对每一列进行变换，最终可以得到图1.1的QR 矩阵。
1.2 程序流程图1.2.1 基于吉文斯变换的QR 分解流程图1.2.2 基于豪斯霍尔德变换的QR 分解流程图计算中间参数计算中间参数()()()()2TTsgn (,):,=,,+1:,R i i R i n i R i i R i n i σωσ=−⎡⎤+⎣⎦T 2,2100u H I uu H H ωω==−⎡⎤=⎢⎥⎣⎦图1.3 基于豪斯霍尔德变换的QR 分解流程图由上述流程图也可以明显的看出来，吉文斯变换的QR 分解需要的循环步骤大于豪斯霍尔德的步骤，与理论分析的一致。
1.3 程序使用说明1.3.1 基于吉文斯变换的QR 分解程序说明基于吉文斯变换的QR 分解程序主要有两个函数构成，主函数如下所示，包含相应的注释，主函数主要就是进行QR 求解，以及方程组的求解。
function [Q,R,x]=Givens(A,b)%基于吉文斯变换的QR分解与方程组求解%输入参数说明：A:参数矩阵% b：列向量m=size(A,1); %矩阵的行数n=size(A,2); %矩阵的列数R=A; %给定QR分解矩阵的初始矩阵RQ=eye(m); %给定QR分解矩阵的初始正交矩阵Qfor i=1:n-1 %分别循环对每一列[a1,a2...an] 求解[P1,P2..Pn] for j=i+1:m %对aij元素消0，求解PijAi=R(:,i); %读取第i列的所有元素Pij=Pij_cal(Ai,i,j); %调用Pij矩阵计算函数Q=Q*Pij'; %求解正交矩阵QR=Pij*R; %求解Rendend%根据QR分解结果求解方程组b=Q'*b; %方程组两边同时左乘正交矩阵Q的逆x(n)=b(n)/R(n,n); %计算初值；for i=n:-1:1 %循环求解x的值sum=0; %每次对求和变量进行初始化for j=i+1:nsum=sum+R(i,j)*x(j);endx(i)=(b(i)-sum)/R(i,i);endend子函数的主要是进行每一个元素的消去的矩阵参数计算，具体如下所示：function [Pij]=Pij_cal(Ai,i,j)%计算每个元素的正交矩阵%输入参数说明：Ai:第i行的列向量% i，j：需要消零元素的编号ai=Ai(i); %需要消元行的首元素赋值；aj=Ai(j); %消元的元素赋值；c=sqrt(ai^2+aj^2); %用于计算角度值；cos_theta=ai/c; %计算cos值sin_theta=aj/c; %计算sin值Pij=eye(length(Ai)); %给定m*m的初始单位矩阵；Pij(i,i)=cos_theta;Pij(i,j)=sin_theta;Pij(j,i)=-sin_theta;Pij(j,j)=cos_theta;end1.3.2 基于豪斯霍尔德变换的QR分解程序说明在进行豪斯霍尔德变换的编程时，我们所需奥解决的主要问题是进行QR矩阵的求解，而对于矩阵的相乘以及矩阵基本的运算和求解向量的范数而言，在进行编程的时候并没有使用自己编写的函数，而是调用了MATLAB中相应的算法。
具体的程序及相应的注释如下所示：function [Q,R,x]=householder(A,b)%基于豪斯霍尔德的QR求解与方程组求解m=size(A,1); %矩阵的行数n=size(A,2); %矩阵的列数R=A; %给定QR分解矩阵的初始矩阵RQ=eye(m); %给定QR分解矩阵的初始正交矩阵Qfor i=1:n-1if R(i,i)>=0 %对zgama符号进行判断sigma=-norm(R(i:n,i));elsesigma=norm(R(i:n,i));endw=[R(i,i)-sigma,R(i+1:n,i)']'; %求解中间计算参数wu=w/norm(w,2);H=eye(m-i+1)-2*u*u'; %求解矩阵HiH=blkdiag(eye(i-1),H); %将矩阵扩展为m*m 的矩阵 Q=Q*H'; %求解Q 矩阵 R=H*R; %求解R 矩阵 end%根据QR 分解结果求解方程组b=Q'*b; %方程组两边同时左乘正交矩阵Q 的逆 x(m)=b(m)/R(m,m); %计算初值；for i=n:-1:1 %循环求解x 的值sum=0; %每次对求和变量进行初始化 for j=i+1:nsum=sum+R(i,j)*x(j); endx(i)=(b(i)-sum)/R(i,i); end end1.4 算例计算结果为了验证所编写的程序是否正确，采用了作业题中第74页的计算实习2.4进行验证，该题具体的数据如式（1.5）所示。
54756753941287886537810987556 , 5791197553688910189587877810105756759101052A b ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦（1.5）将相应的数据输入到MATLAB 中，分别运行基于吉文斯变换的QR 分解程序“[Q,R,x]=Givens(A,b)”和豪斯霍尔德程序“[Q,R,x]=householder(A,b)”可以得到相应的方程组的解如式（1.6）所示：[]T1111111x = （1.6）经过比较发现，采用两种算法得到的QR 分解矩阵，如式（1.7）所示，经过比较发现两者得到的结果在正负号方面存在差异，主要是豪斯霍尔德算法中对于σ取值确定上存在正负的情况，因此略有不同。
0.33330.33720.28770.43670.09820.22220.66530.26670.91200.01360.24620.11810.07060.13190.46670.14980.19610.14070.57040.44460.42250.33330.06570.56920.60350.05840.42360.11660.40000.0251givensQ −−−−−−−−−−−=−0.13600.37780.38150.62130.38290.46670.14980.72940.15600.15470.42000.05500.33330.06860.06400.44450.69150.07780.44520.33330.33720.28770.43670.09820.22220.6householderQ ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥−−⎢⎥−−−−⎢⎥⎢⎥−−−⎣⎦−−−=6530.26670.91200.01360.24620.11810.07060.13190.46670.14980.19610.14070.57040.44460.42250.33330.06570.56920.60350.05840.42360.11660.40000.02510.13600.37780.38150.62130.38290.46670.149−−−−−−−−−−−−−−−−−−80.72940.15600.15470.42000.05500.33330.06860.06400.44450.69150.07780.4452⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥−⎢⎥⎢⎥−−−⎣⎦得到的上三角R 矩阵的结果如下所示15.000019.533320.933319.933321.600021.266719.80000.00007.4464 2.6996 2.90553.0995 2.3624 1.10660.00000.0000 3.2416 3.3580 1.68830.4811 2.30360.00000.00000.0000 3.73420.7405 1.6506 1.11390.00000.givensR −−−−−=−−−−householder00000.00000.0000 3.2303 3.2048 3.90190.00000.00000.00000.00000.0000 1.98010.07380.00000.00000.00000.00000.00000.00000.978615.000019.533320.933319.93321.60R ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥−−⎢⎥−−−−⎢⎥⎢⎥−−−−⎣⎦−−−−−=0021.266719.80000.00007.4464 2.6996 2.9055 3.0995 2.3624 1.10660.00000.0000 3.2416 3.3580 1.68830.4811 2.30360.00000.00000.0000 3.73420.7405 1.6506 1.11390.00000.00000.00000.0000 3.2303 3.2048 3.−−−−−−−−−−−−−−−−−−−−90190.00000.00000.00000.00000.0000 1.98010.07380.00000.00000.00000.00000.00000.00000.9786⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥−−−⎢⎥⎢⎥⎣⎦为了验证所计算结果是否正确，采用MATLAB 中关于QR 分解的库算法，在命令栏中输入“[Q,R]=qr(A)”得到的结果与采用豪斯霍尔德的计算结果相同，在命令栏中输入“[x]=inv(A)*b ”得到的方程组的解与计算结果相同，验证了算法的正确性。