资源描述
word工业大学计算机与信息学院实验报告课 程：计算方法专业班级： 学 号： 姓 名： 33 / 33Java界面其实都不难按照程序流程图就可以完成了实验一插值与拟合一、 实验目的(1) 明确插值多项式和分段插值多项式各自的优缺点；(2) 编程实现三次样条插值算法，分析实验结果体会高次插值产生的龙格现象；(3) 理解最小二乘拟合，并编程实现线性拟合，掌握非线性拟合转化为线性拟合的方法(4) 运用常用的插值和拟合方法解决实际问题。二、 实验容(1)对于f(x)=1/(1+x*x)实现三次样条插值(2)实现最小二乘法的直线拟合数据如下：165123150123141187126172125148三、 根本原理计算公式(1)三次样条插值在每个节点上具有2阶导数。(2) 最小二乘法拟合直线为y=a+bx，而a，b有如下等式N为给出的数据点的总个数 ； 四、算法设计与实现流程图，关键点最小二乘法直线拟合：输入数据后，按照公式计算a，b。用得到的拟合直线计算预测点的近似函数值。五、输入与输出(1)三次样条插值输入：区间长度，n+1个数据点，预测点输出：预测点的近似函数值，准确值，与误差(2)最小二乘法直线拟合输入：n个数据点，预测点输出：预测点的近似函数值六、结果讨论和分析代码三次样条插值#include#include #define N 10using namespace std;double u0(double x)return (x-1)*(x-1)*(2*x+1);double u1(double x)return x*x*(3-2*x);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; fin.open(E:t.txt); if(!fin) couterror opening input streamendl; system(pause); return 0; double xN+1,yN+1,mN+1,AN,BN,CN;double hN;double aN,bN;double f0,fn;double temp;int i;for(i=0;ixiyi;finf0fn;h0=x1-x0;for(i=1;iN;i+)hi=xi+1-xi;ai=hi-1/(hi-1+hi);bi=3*(1-ai)*(yi-yi-1)/hi-1+ai*(yi+1-yi)/hi);m1=b1-(1-a1)*f0;mN-1=bN-1-aN-1*fn;for(i=2;iN-1;i+)mi=bi; for(i=1;iN;i+)Bi=2;Ci=ai;for(i=2;i0;i-)mi=mi-Ci*mi+1;coutplease:(输入插值节点在x0到xN围)temp) double tt=temp;if(tempxN)cout插值节点为tt超出插值围endl;continue;for(i=1;i=N;i+)if(tempxi)break;temp=(temp-xi-1)/hi-1; temp=s3(temp,yi-1,yi,mi-1,mi,hi-1); cout插值节点为tt准确值为f(tt)插值结果为temp误差为f(tt)-tempendl; system(pause); fin.close();return 0;最小二乘法的直线拟合#include#include#define n 5using namespace std;double sum(double x,int k) int i; double sum=0; for(i=0;ik;i+) sum=sum+xi; return sum; double sum2(double x,int k) int i; double sum=0; for(i=0;ik;i+) sum=sum+xi*xi; return sum; double sumxy(double x,double y,int k) int i; double sum=0; for(i=0;ik;i+) sum=sum+xi*yi; return sum; int main() ifstream fin; fin.open(E:t.txt); if(!fin) couterror opening input streamendl; system(pause); return 0; double xn,yn,a,b; double x0,y0; int i; for(i=0;ixiyi; 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: a,b: b,拟合直线为y=a+bxendl; coutx0) y0=a+b*x0; cout当x=x0,y=y0endl; cout请输入插值节点x：; system(pause); fin.close(); return 0;实验二数值积分一、 实验目的(1) 熟悉复化梯形方法、复化Simpson方法、梯形递推算法、龙贝格算法；(2) 能编程实现龙贝格算法和中点加速；(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四、 算法设计与实现流程图，关键点中点加速梯形递推算法流程图龙贝格算法流程图：输入数据后根据公式计算导数值五、输入与输出(1) 用龙贝格算法计算输入：积分区间，误差限输出：序列Tn，Sn,Rn与积分结果(2) 用中点加速方法计算的一阶导数输入：求导节点，步长 输出：求得的导数值，准确值六、结果讨论和分析代码龙贝格算法#include#include#includeusing namespace std;double f(double x) if(x=0)return 1; return sin(x)/x; int main() ifstream fin; fin.open(E:t.txt); if(!fin) couterror opening input streamabe; cout积分区间为a,b,要求精度为eendl; coutk T2 S2 C2 R2endl; h=b-a; t1=(f(a)+f(b)*h/2; cout0 t1endl; int k; for(k=1;k=10;k+,h=h/2,t1=t2,s1=s2) s=0; x=a+h/2; do s=s+f(x); x=x+h; while(xb); t2=t1/2+h*s/2; s2=t2+(t2-t1)/3; if(k=1) coutk t2 s2endl; continue; c2=s2+(s2-s1)/15; if(k=2) coutk t2 s2 c2endl; c1=c2; continue; r2=c2+(c2-c1)/63; coutk t2 s2 c2r2endl; if(k=3) r1=r2; c1=c2; continue; if(fabs(r2-r1)e) cout数值积分结果为r2endl; break; r1=r2; c1=c2; system(pause); return 0;中点加速算法#include#include#includeusing namespace std;double f(double x) return exp(x); double f1(double 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,double h) return 64*g2(x,h/2)/63-g2(x,h)/63; int main() ifstream fin; fin.open(E:t.txt); if(!fin) couterror opening input streamah) cout当x=a ,步长h= h,x处一阶导数值准确值为f1(a),中点加速求得x处一阶导数值为g3(a,h),误差为f1(a)-g3(a,h)endl; system(pause); fin.close(); return 0;实验三非线性方程求根迭代法一、 实验目的(1) 熟悉非线性方程求根简单迭代法，牛顿迭代与牛顿下山法(2) 能编程实现牛顿下山法(3) 认识选择迭代格式的重要性(4) 对迭代速度建立感性的认识；分析实验结果体会初值对迭代的影响二、 实验容用牛顿下山法解方程初值为0.6三、 根本原理计算公式求非线性方程组的解是科学计算常遇到的问题，有很多实际背景各种算法层出不穷，其中迭代是主流算法。只有建立有效的迭代格式，迭代数列才可以收敛于所求的根。因此设计算法之前，对于一般迭代进展收敛性的判断是至关重要的。牛顿法也叫切线法，是迭代算法中典型方法，只要初值选取适当，在单根附近，牛顿法收敛速度很快，初值对于牛顿迭代至关重要。当初值选取不当可以采用牛顿下山算法进展纠正。一般迭代：牛顿公式：牛顿下山公式：下山因子下山条件四、算法设计与实现流程图，关键点牛顿下山算法流程图五、输入与输出输入：初值，误差限，迭代最大次数，下山最大次数输出：近似根各步下山因子 六、结果讨论和分析代码牛顿下山法#include#include#includeusing namespace std;double f(double x)/函数 return x*x*x-x-1; double f1(double x)/一阶导数 return 3*x*x-1; int main() ifstream fin; fin.open(E:t.txt ); if(!fin) cout Error opening input stream x0eMN; cout迭代初值为x0,误差限为e,最大下山次数为M,最大迭代次数为Nendl; x1=0; for(k=0;kN&f1(x0)!=0;k+) i=0; j=1; temp=f(x0)/f1(x0); do x1=x0-temp/j; coutx0:x0,x1:x1,下山次数:i,下山因子:1/j,迭代次数:k=M)break; while(fabs(f(x1)=fabs(f(x0);if(i=M) cout重新选择x0endl; break; if(fabs(x1-x0)e) cout迭代成功,求得结果为x1=N)cout迭代失败endl; if(f1(x0)=0)cout一阶导数为0endl; fin.close(); system(pause); return 0;实验四求解线性方程组一、 实验目的(1) 熟悉求解线性方程组的有关理论和方法；(2) 能编程实现高斯-塞德尔迭代法、列主元高斯消去法、LU分解法(3) 通过测试，进一步了解各种方法的优缺点(4) 根据不同类型的方程组，选择适宜的数值方法二、 实验容(1) 用高斯-塞德尔迭代法求(2) 列主元高斯消去法求(3) LU分解法求解方程组三、 根本原理计算公式线性方程组大致分迭代法和直接法。只有收敛条件满足时，才可以进展迭代。雅可比与高斯-塞德尔是最根本的两类迭代方法，最大区别是迭代过程中是否引用新值进展剩下的计算。消元是最简单的直接法，并且也十分有效的，列主元高斯消去法对求解一般的线性方程组都适用，同时可以用来求矩阵对应的行列式。约当消去实质是经过初等行变换将系数矩阵化为单位阵，主要用来求矩阵的逆。在使用直接法，要注意从空间、时间两方面对算法进展优化。高斯-塞德尔迭代：列主元高斯消去法：列主元消元 回代四、算法设计与实现流程图，关键点列主元的高斯消去流程图G-S迭代算法流程图LU分解法：依次求得L、U、y和x五、输入与输出(1) 用高斯-塞德尔迭代法输入：系数矩阵A，最大迭代次数N，初始向量，误差限e输出：解向量(2) 列主元高斯消去法输入：系数矩阵A输出：解向量(3) LU分解法输入：系数矩阵A输出：解向量六、结果讨论和分析代码高斯塞德尔迭代#include#include#include#define n 3using namespace std;void show(double an+1n+1,double bn+1) int i,j; cout原方程为：endl; for (i = 1; i = n; i+) for (j = 1; j n; j+)if (aij + 1 0)coutaij*xj;elsecoutaij*x j+;coutaij*x j=biendl; int main() ifstream fin; fin.open(E:t.txt ); if(!fin) cout Error opening input stream eN; for(i=1;ixi; yi=xi; cout初始向量为:; for(i=1;in;i+)coutxi,; coutxiendl; for(i=1;i=n;i+) for(j=1;jaij; finbi; 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+aij*yj; yi=(bi-temp)/aii; max=fabs(y1-x1); for(i=2;i=n;i+) if(maxfabs(yi-xi) max=fabs(yi-xi); if(maxe) cout高斯赛德尔迭代求得原方程的解为endl; for(i=1;i=n;i+)coutxi为yi ; coutendl; break; if(k=N) cout迭代失败endl; break; else k+; for(i=1;i=n;i+)xi=yi; fin.close(); system(pause); return 0;高斯消去#include#include#include#define n 3using namespace std;void show(double an+1n+1,double bn+1) int i,j; cout原方程为：endl; for (i = 1; i = n; i+) for (j = 1; j n; j+)if (aij + 1 0)coutaij*xj;elsecoutaij*x j+;coutaij*x j=biendl; int main() ifstream fin; fin.open(E:t.txt); if(!fin) couterror opening input streamendl; system(pause); return 0; double an+1n+1,bn+1,d,t; int i,j,k,l; for(i=1;i=n;i+) for(j=1;jaij; finbi; show(a,b); k=1; do d=akk; l=k; i=k+1; do if(in)break; if(fabs(aik)fabs(d) d=aik; l=i; if(i=n)break; i+; while(true); if(d=0) cout奇异endl; system(pause); fin.close(); return 0; if(l!=k) for(j=k;j=n;j+) t=alj; alj=akj; akj=t; t=bk; bk=bl; bl=t; for(j=k+1;j=n;j+) akj=akj/akk; bk=bk/akk; for(i=k+1;i=n;i+) for(j=k+1;j=n;j+) aij=aij-aik*akj; for(i=k+1;i=1;i-) t=0; for(j=i+1;j=n;j+) t=t+aij*bj; bi=bi-t; cout列主元的高斯消去法求得原方程的解为: ; for(i=1;i=n;i+)coutxi为bi ; coutendl; system(pause); fin.close(); return 0;LU分解#include#include#define n 3using namespace std;void show(double an+1n+1,double bn+1) int i,j; cout原方程为：endl; for (i = 1; i = n; i+) for (j = 1; j n; j+)if (aij + 1 0)coutaij*xj;elsecoutaij*x j+;coutaij*x j=biendl; int main() ifstream fin; fin.open(E:t.txt); if(!fin) couterror opening input streamendl; system(pause); return 0; double an+1n+1,bn+1,un+1n+1,ln+1n+1,xn+1,yn+1; double t; int i,j,k; for(i=1;i=n;i+) for(j=1;jaij; finbi; show(a,b); for(i=1;i=n;i+) for(j=1;ji;j+) t=0; for(k=1;kj;k+) t=t+lik*ukj; lij=(aij-t)/ujj; for(j=i;j=n;j+) t=0; for(k=1;ki;k+) t=t+lik*ukj; uij=aij-t; 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+uij*xj; xi=(yi-t)/uii; coutLU分解法求得原方程的解为: ; for(i=1;i=n;i+)coutxi为xi ; coutendl; system(pause); fin.close(); return 0;实验五数值微分一、 实验目的(1) 熟悉数值微分中Euler法，改良Euler法，Rung-Kutta方法；(2) 能编程实现亚当姆斯方法，Rung-Kutta方法；(3) 通过实验结果分析各个算法的优缺点;(4) 明确步长对算法的影响并理解变步长的Rung-Kutta方法二、 实验容0 x1取h=0.1时用亚当姆斯方法，Rung-Kutta方法求其数值解并与准确解进展比拟。三、 根本原理计算公式在许多科学技术问题中，建立的模型常常以常微分方程的形式表示。然而，除了少数特殊类型的常微分方程能用解析方法求其准确解外，要给出一般方程解析解的表达式是困难的。所以只能用近似方法求其数值解，在实际工作中常用计算机求常微分方程的数值解。所谓常微分方程的数值解即对于常微分方程初值问题计算出在一系列节点 a = x0 x1 xn= b 处的未知函数 y(x)近似值y0,y1,yn，即找到一系列离散点x0,y0x1,y1x2,y2xn,yn近似满足常微分方程。数值解法的根本思想用差商代替导数,实现连续问题离散化，选取不同的差商代替导数可以得到不同公式。改良欧拉公式是常用方法之一，包括预测和校正两步。先用欧拉公式进展预报，再将预报值代入梯形公式进展校正，从而达到二阶精度。通过龙格-库塔法我们可以获得更高精度。经典龙格-库塔法即在区间xn,xn+1取四点，并对这四点的斜率进展加权平均作为平均斜率，通过泰勒公式寻找使局部截断误差为O(h5)即4阶精度的参数满足条件。四阶经典龙格-库塔公式四、 算法设计与实现流程图，关键点亚当姆斯方法经典龙格库塔算法五、 输入与输出输入：求解区间，初值，数值解个数输出：数值解六、 结果讨论和分析代码龙格-库塔方法#include#include#includeusing namespace std;double f(double x,double y) return y-2*x/y;/函数 double f0(double x) return pow(1+2*x,0.5); int main() ifstream fin; fin.open(E:t.txt); if(!fin) couterror opening input streamx0y0hN; cout龙格-库塔方法求解常微分方程,初始条件为步长h= h ,初值为x0= x0 ,y0= y0endl; n=1; do x1=x0+h; k1=f(x0,y0); k2=f(x0+h/2,y0+h*k1/2); k3=f(x0+h/2,y0+h*k2/2); k4=f(x1,y0+h*k3); y1=y0+h*(k1+2*k2+2*k3+k4)/6; cout当x= x1 时,y准确值为f0(x1),用龙格-库塔求得y为y1,误差为f0(x1)-y1endl; n+; x0=x1; y0=y1; while(n=N); system(pause); fin.close(); return 0;亚当姆斯方法#include#include#includeusing namespace 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;in-1;i+) k1=f(xi,yi); k2=f(xi+h/2,yi+h*k1/2); k3=f(xi+h/2,yi+h*k2/2); k4=f(xi+h,yi+h*k3); yi+1=yi+h*(k1+2*k2+2*k3+k4)/6; int main() ifstream fin; fin.open(E:t.txt); if(!fin) couterror opening input streamx0y0hN; cout亚当姆斯方法求解常微分方程,初始条件为步长h= h ,初值为x0= x0 ,y0= y0endl; for(k=0;k4;k+)xk=x0+k*h; fc(y,x,4,h); for(k=0;k4;k+)y1k=f(xk,yk); for(k=1;k4;k+)cout当x= xk 时,y准确值为f0(xk),用亚当姆斯求得y为yk,误差为f0(xk)-ykendl; n=4; do x4=x3+h; yp=y3+h*(55*y13-59*y12+37*y11-9*y10)/24; yp1=f(x4,yp); y4=y3+h*(9*yp1+19*y13-5*y12+y11)/24; y14=f(x4,y4); cout当x= x4 时,y准确值为f0(x4),用亚当姆斯求得y为y4,误差为f0(x4)-y4endl; n+; x3=x4; y3=y4; for(k=0;k4;k+)y1k=y1k+1; while(n=N); system(pause); fin.close(); return 0;
展开阅读全文