合肥工业大学 计算机专业 计算方法实验报告材料

《合肥工业大学 计算机专业 计算方法实验报告材料》由会员分享，可在线阅读，更多相关《合肥工业大学 计算机专业 计算方法实验报告材料（33页珍藏版）》请在装配图网上搜索。

1、word 工业大学 计算机与信息学院 实验报告 课 程：计算方法 专业班级： 学 号： 姓 名： 33 / 33 Java界面 其实都不难按照程序流程图就可以完成了 实验一插值与拟合 一、 实验目的 (1) 明确插值多项式和分段插值多项式各自的优缺点； (2) 编程实现三次样条插值算法，分析实验结果体会高次插值产生的龙格现象； (3) 理解最小二乘拟合，并编程实现线性拟合，掌握非线性拟合转化为线性拟合的方法 (4) 运用常用的插值和拟合方法

2、解决实际问题。 二、 实验容 (1)对于f(x)=1/(1+x*x)实现三次样条插值 (2)实现最小二乘法的直线拟合 数据如下： 165 123 150 123 141 187 126 172 125 148 三、 根本原理〔计算公式〕 (1)三次样条插值在每个节点上具有2阶导数。 (2) 最小二乘法拟合直线为y=a+bx，而a，b有如下等式〔N为给出的数据点的总个数〕 ； 四、算法设计与实现〔流程图，关键点〕 最小二乘法直线拟合：输入数据后，按照公式计算a，b。用得到的拟合直线计算预测点的近似函数值。 五、输入与输出

3、(1)三次样条插值 输入：区间长度，n+1个数据点，预测点 输出：预测点的近似函数值，准确值，与误差 (2)最小二乘法直线拟合 输入：n个数据点，预测点 输出：预测点的近似函数值 六、结果讨论和分析 代码 三次样条插值 #include #include #define N 10 using namespace std; double u0(double x){ return (x-1)*(x-1)*(2*x+1); } double u1(double x){ return x*x*(3-2*x); }

4、 double v0(double x){ return x*(x-1)*(x-1); } double v1(double x){ return x*x*(x-1); } double s3(double x,double y,double y1,double m,double m1,double h){ return u0(x)*y+u1(x)*y1+h*v0(x)*m+h*v1(x)*m1; } double f(double x){ return 1/(1+x*x); } int main(){ ifstream fin;

5、 fin.open("E:\\t.txt"); if(!fin){ cout<<"error opening input stream"<

6、N;i++){ fin>>x[i]>>y[i]; } fin>>f0>>fn; h[0]=x[1]-x[0]; for(i=1;i

7、){ B[i]=2; C[i]=a[i]; } for(i=2;i0;i--){ m[i]=m[i]-C[i]*m[i+1]; }

8、cout<<"please:(输入插值节点在"<>temp){ double tt=temp; if(tempx[N]){ cout<<"插值节点为"<

9、i],h[i-1]); cout<<"插值节点为"< #include #define n 5 using namespace std; double sum(double x[],int k){ int i; double sum=

10、0; for(i=0;i

11、uble sum=0; for(i=0;i

12、 } double x[n],y[n],a,b; double x0,y0; int i; for(i=0;i>x[i]>>y[i]; } b=(n*sumxy(x,y,n)-sum(x,n)*sum(y,n))/(n*sum2(x,n)-sum(x,n)*sum(x,n)); a=(sum(y,n)-b*sum(x,n))/n; cout<<"最小二乘法直线拟合得到a: "<

13、合直线为y="<>x0){ y0=a+b*x0; cout<<"当x="<

14、点加速； (3) 理解并掌握自适应算法和收敛加速算法的根本思想； (4) 分析实验结果体会各种方法的准确度，建立计算机求解定积分问题的感性认识 二、 实验容 (1) 用龙贝格算法计算 (2) 用中点加速方法计算的一阶导数 三、 根本原理〔计算公式〕 (1)龙贝格算法 梯形递推公式 加权平均公式： (2) 中点加速 中点公式: G(h)=(f(a+h)-f(a-h))/2/h 加权平均:G1(h)=4*G(h/2)/3-G(h)/3 G2(h)=16*G1(h/2)/15-G1(h)/15 G3(h)=64*G2(h/2)/63-G2(h)/63 四

15、、 算法设计与实现〔流程图，关键点〕 中点加速梯形递推算法流程图 龙贝格算法流程图 ：输入数据后根据公式计算导数值 五、输入与输出 (1) 用龙贝格算法计算 输入：积分区间，误差限 输出：序列Tn，Sn,,Rn与积分结果 (2) 用中点加速方法计算的一阶导数 输入：求导节点，步长 输出：求得的导数值，准确值 六、结果讨论和分析 代码 龙贝格算法 #include #include #include using namespace std; double f(double x){

16、 if(x==0)return 1; return sin(x)/x; } int main(){ ifstream fin; fin.open("E:\\t.txt"); if(!fin){ cout<<"error opening input stream"<

17、,r2; double x,h,s; fin>>a>>b>>e; cout<<"积分区间为["<

18、o{ s=s+f(x); x=x+h; }while(x

19、c2<

20、"数值积分结果为"< #include #include using namespace std; double f(double x){ return exp(x); } double f1(doub

21、le x){ return exp(x); } double g(double x,double h){ return (f(x+h)-f(x-h))/2/h; } double g1(double x,double h){ return 4*g(x,h/2)/3-g(x,h)/3; } double g2(double x,double h){ return 16*g1(x,h/2)/15-g1(x,h)/15; } double g3(double x,d

22、ouble h){ return 64*g2(x,h/2)/63-g2(x,h)/63; } int main(){ ifstream fin; fin.open("E:\\t.txt"); if(!fin){ cout<<"error opening input stream"<>a>

23、>h) cout<<"当x="<

24、结果体会初值对迭代的影响 二、 实验容 用牛顿下山法解方程〔初值为0.6〕 三、 根本原理〔计算公式〕 求非线性方程组的解是科学计算常遇到的问题，有很多实际背景．各种算法层出不穷，其中迭代是主流算法。只有建立有效的迭代格式，迭代数列才可以收敛于所求的根。因此设计算法之前，对于一般迭代进展收敛性的判断是至关重要的。牛顿法也叫切线法，是迭代算法中典型方法，只要初值选取适当，在单根附近，牛顿法收敛速度很快，初值对于牛顿迭代 至关重要。当初值选取不当可以采用牛顿下山算法进展纠正。 一般迭代： 牛顿公式： 牛顿下山公式： 下山因子 下山条件 四、算法设计与实现〔流程图，

25、关键点〕 牛顿下山算法流程图 五、输入与输出 输入：初值，误差限，迭代最大次数，下山最大次数 输出：近似根各步下山因子 六、结果讨论和分析 代码 牛顿下山法 #include #include #include using namespace std; double f(double x){//函数 return x*x*x-x-1; } double f1(double x){//一阶导数 return 3*x*x-1; }

26、 int main(){ ifstream fin; fin.open("E:\\t.txt" ); if(!fin){ cout << "Error opening input stream" << endl; system("pause"); return 0; } double x0,x1,e,j,temp; int N,M,k,i; fin>>x0>>e>>M>>N; cout<<"迭代初值为"<

27、e<<",最大下山次数为"<=M)break;

28、 }while(fabs(f(x1))>=fabs(f(x0))); if(i>=M){ cout<<"重新选择x0"<

29、 x0=x1; } if(k>=N)cout<<"迭代失败"<

30、尔迭代法求 (2) 列主元高斯消去法求 (3) LU分解法求解方程组 三、 根本原理〔计算公式〕 线性方程组大致分迭代法和直接法。只有收敛条件满足时，才可以进展迭代。雅可比与高斯-塞德尔是最根本的两类迭代方法，最大区别是迭代过程中是否引用新值进展剩下的计算。消元是最简单的直接法，并且也十分有效的，列主元高斯消去法对求解一般的线性方程组都适用，同时可以用来求矩阵对应的行列式。约当消去实质是经过初等行变换将系数矩阵化为单位阵，主要用来求矩阵的逆。在使用直接法，要注意从空间、时间两方面对算法进展优化。 高斯-塞德尔迭代： 列主元高斯消去法：列主元 消元 回代

31、 四、算法设计与实现〔流程图，关键点〕 列主元的高斯消去流程图 G-S迭代算法流程图 LU分解法：依次求得L、U、y和x 五、输入与输出 (1) 用高斯-塞德尔迭代法 输入：系数矩阵A，最大迭代次数N，初始向量，误差限e 输出：解向量 (2) 列主元高斯消去法 输入：系数矩阵A 输出：解向量 (3) LU分解法 输入：系数矩阵A 输出：解向量 六、结果讨论和分析 代码 高斯塞德尔迭代 #include #include #

32、include #define n 3 using namespace std; void show(double a[n+1][n+1],double b[n+1]){ int i,j; cout<<"原方程为："<

33、<"*x"<< j<<"="<

34、,N; fin>>e>>N; for(i=1;i<=n;i++){ fin>>x[i]; y[i]=x[i]; } cout<<"初始向量为:"; for(i=1;i>a[i

35、][j]; fin>>b[i]; } show(a,b); k=0; while(true){ for(i=1;i<=n;i++){ temp=0; for(j=1;j<=n;j++) if(j!=i)temp=temp+a[i][j]*y[j]; y[i]=(b[i]-temp)/a[i][i]; } max=fabs(y[1]-x[1]); for(i=2;i<=n;i+

36、+) if(max

37、 k++; for(i=1;i<=n;i++)x[i]=y[i]; } } fin.close(); system("pause"); return 0; } 高斯消去 #include #include #include #define n 3 using namespace std; void show(double a[n+1][n+1],double b[n+1]){ int i,j; cout<<"原方程为："<

38、dl; for (i = 1; i <= n; i++) { for (j = 1; j < n; j++) if (a[i][j + 1] < 0) cout<

39、"error opening input stream"<>a[i][j]; fin>>b[i]; }

40、show(a,b); k=1; do{ d=a[k][k]; l=k; i=k+1; do{ if(i>n)break; if(fabs(a[i][k])>fabs(d)){ d=a[i][k];

41、 l=i; } if(i==n)break; i++; }while(true); if(d==0){ cout<<"奇异"<

42、 system("pause"); fin.close(); return 0; } if(l!=k){ for(j=k;j<=n;j++){ t=a[l

43、][j]; a[l][j]=a[k][j]; a[k][j]=t; } t=b[k]; b[k]=b[l];

44、 b[l]=t; } for(j=k+1;j<=n;j++){ a[k][j]=a[k][j]/a[k][k]; } b[k]=b[k]/a[k][k]; for(i=k+1;i<=n;i++){

45、 for(j=k+1;j<=n;j++){ a[i][j]=a[i][j]-a[i][k]*a[k][j]; } } for(i=k+1;i<=n;i++)b[i]=b[i]-a[i

46、][k]*b[k]; if(k==n)break; k++; }while(true); for(i=n-1;i>=1;i--){ t=0; for(j=i+1;j<=n;j++){ t=t+a[i][j]*b[j];

47、 } b[i]=b[i]-t; } cout<<"列主元的高斯消去法求得原方程的解为:" ; for(i=1;i<=n;i++)cout<<"x"<

48、eam> #include #define n 3 using namespace std; void show(double a[n+1][n+1],double b[n+1]){ int i,j; cout<<"原方程为："<

49、

50、][n+1],x[n+1],y[n+1]; double t; int i,j,k; for(i=1;i<=n;i++){ for(j=1;j<=n;j++)fin>>a[i][j]; fin>>b[i]; } show(a,b); for(i=1;i<=n;i++){ for(j=1;j

51、 t=0; for(k=1;k

52、=i;j<=n;j++){ t=0; for(k=1;k

53、 } for(i=1;i<=n;i++){ t=0; for(j=1;j=1;i--){ t=0; for(j=i+1;j<=n;j++)t=t+u[i][j]*x[j];

54、 x[i]=(y[i]-t)/u[i][i]; } cout<<"LU分解法求得原方程的解为:" ; for(i=1;i<=n;i++)cout<<"x"<

55、ung-Kutta方法； (3) 通过实验结果分析各个算法的优缺点; (4) 明确步长对算法的影响并理解变步长的Rung-Kutta方法 二、 实验容 0 < x<1 取h=0.1时用亚当姆斯方法，Rung-Kutta方法求其数值解并与准确解进展比拟。 三、 根本原理〔计算公式〕 在许多科学技术问题中，建立的模型常常以常微分方程的形式表示。然而，除了少数特殊类型的常微分方程能用解析方法求其准确解外，要给出一般方程解析解的表达式是困难的。所以只能用近似方法求其数值解，在实际工作中常用计算机求常微分方程的数值解。所谓常微分方程的数值解即对于常微分方程初值问题计算出在一系列节点

56、 a = x0< x1<…< xn= b 处的未知函数 y(x)近似值y0,y1,…yn，即找到一系列离散点〔x0,y0〕〔x1,y1〕〔x2,y2〕…〔xn,yn〕近似满足常微分方程。数值解法的根本思想用差商代替导数,实现连续问题离散化，选取不同的差商代替导数可以得到不同公式。 改良欧拉公式是常用方法之一，包括预测和校正两步。先用欧拉公式进展预报，再将预报值代入梯形公式进展校正，从而达到二阶精度。通过龙格-库塔法我们可以获得更高精度。经典龙格-库塔法即在区间[xn,xn+1]取四点，并对这四点的斜率进展加权平均作为平均斜率，通过泰勒公式寻找使局部截断误差为O(h5)〔即4阶精度〕的参数满足

57、条件。 四阶〔经典〕龙格-库塔公式 四、 算法设计与实现〔流程图，关键点〕 亚当姆斯方法 经典龙格库塔算法 五、 输入与输出 输入：求解区间，初值，数值解个数 输出：数值解 六、 结果讨论和分析 代码 龙格-库塔方法 #include #include #include using namespace std; double f(double x,double y){ return y-2*x/y;//函数 } double f

58、0(double x){ return pow(1+2*x,0.5); } int main(){ ifstream fin; fin.open("E:\\t.txt"); if(!fin){ cout<<"error opening input stream"<

59、int N,n; fin>>x0>>y0>>h>>N; cout<<"龙格-库塔方法求解常微分方程,初始条件为步长h="<< h << ",初值为x0=" < #include #include using namesp

61、ace std; double f(double x,double y){ return y-2*x/y;//函数 } double f0(double x){ return pow(1+2*x,0.5); } void fc(double y[],double x[],int n,double h){ double k1,k2,k3,k4; int i; for(i=0;i

62、x[i],y[i]); k2=f(x[i]+h/2,y[i]+h*k1/2); k3=f(x[i]+h/2,y[i]+h*k2/2); k4=f(x[i]+h,y[i]+h*k3); y[i+1]=y[i]+h*(k1+2*k2+2*k3+k4)/6; } } int main(){ ifstream

63、 fin; fin.open("E:\\t.txt"); if(!fin){ cout<<"error opening input stream"<>x[0]>>y[0]>>h>>N; cout<<"亚当姆斯方法求解常微分方程,初始条件为步长h="<<

64、 h << ",初值为x0=" <

65、 yp=y[3]+h*(55*y1[3]-59*y1[2]+37*y1[1]-9*y1[0])/24; yp1=f(x[4],yp); y[4]=y[3]+h*(9*yp1+19*y1[3]-5*y1[2]+y1[1])/24; y1[4]=f(x[4],y[4]); cout<<"当x= "<