计算方法实验报告(附代码).doc

实验一 牛顿下山法实验说明：求非线性方程组的解是科学计算常遇到的问题，有很多实际背景各种算法层出不穷，其中迭代是主流算法。只有建立有效的迭代格式，迭代数列才可以收敛于所求的根。因此设计算法之前，对于一般迭代进行收敛性的判断是至关重要的。牛顿法也叫切线法，是迭代算法中典型方法，只要初值选取适当，在单根附近，牛顿法收敛速度很快，初值对于牛顿迭代至关重要。当初值选取不当可以采用牛顿下山算法进行纠正。牛顿下山公式： 下山因子下山条件实验代码：#include#include#includeusing namespace std;double newton_downhill(double x0,double x1); /牛顿下山法函数，返回下山成功后的修正初值 double Y; /定义下山因子Ydouble k; /k为下山因子Y允许的最小值double dfun(double x)return 3*x*x-1; /dfun()计算f(x)的导数值double fun1(double x)return x*x*x-x-1; /fun1()计算f(x)的函数值double fun2(double x) return x-fun1(x)/dfun(x); /fun2()计算迭代值int N; /N记录迭代次数double e; /e表示要求的精度int main()double x0,x1;coutx0;coute;N=1;if(dfun(x0)=0)coutf(x0)=0,无法进行牛顿迭代!endl;x1=fun2(x0);coutx0setw(18)x1setw(18)esetw(25)f(x1)-f(x0)endl; coutsetiosflags(ios:fixed)setprecision(6)x0 x1 fabs(x1-x0) fabs(fun1(x1)-fabs(fun1(x0)=fabs(fun1(x0) /初值不满足要求时,转入牛顿下山法x1=newton_downhill(x0,x1); /牛顿下山法结束后,转入牛顿迭代法进行计算while(fabs(x1-x0)=e) /当精度不满足要求时N=N+1;x0=x1; if(dfun(x0)=0)cout迭代途中f(x0)=0,无法进行牛顿迭代!endl; x1=fun2(x0);coutsetiosflags(ios:fixed)setprecision(6)x0 x1 fabs(x1-x0)endl;cout迭代值为:setiosflags(ios:fixed)setprecision(6)x1n;cout迭代次数为:Nendl;return 0;double newton_downhill(double x0,double x1)Y=1;coutk;while(fabs(fun1(x1)=fabs(fun1(x0)if(Yk)Y=Y/2;else cout下山失败!;exit(0);x1=x0-Y*fun1(x0)/dfun(x0);/下山成功则cout下山成功!Y=Y,转入牛顿迭代法计算!endl;return x1;实验结果：实验二 高斯-塞德尔迭代法实验说明：线性方程组大致分迭代法和直接法。只有收敛条件满足时，才可以进行迭代。雅可比及高斯-塞德尔是最基本的两类迭代方法，最大区别是迭代过程中是否引用新值进行剩下的计算。图4.1G-S迭代算法流程图高斯-塞德尔迭代： 实验代码： #include #include #include using namespace std; class gsdl private: int n; double *a, *b, *x, eps; public: gsdl (int nn) int i; n = nn; a = new double*n; /动态分配内存空间 for (i=0; in; i+) ai = new doublen; b = new doublen; x = new doublen; void input (); /从文件读入系数矩阵A、常数向量B以及e void a_gsdl (); void output (); /输出结果到文件并显示 gsdl () int i; for (i=0; in; i+) delete ai; delete a; delete b, x; ; void gsdl:input () /从文件读入系数矩阵A、常数向量B以及e int i, j; char str120; cout str1; ifstream fin (str1); if (!fin) cout n不能打开这个文件 str1 endl; exit(1); for (i=0; in; i+) /读入矩阵A for (j=0; jaij; for (i=0; ibi; /读入常数向量B fin eps; /读入e fin.close (); void gsdl:a_gsdl () int i,j; double p,t,s,q; for (i=0; i=n-1; i+) p=0.0; xi=0.0; for (j=0; j=fabs(aii) cout n程序工作失败！ =eps) p=0.0; for (i=0; i=n-1; i+) t=xi; s=0.0; for (j=0; jp) p=q; void gsdl:output () /输出结果到文件并显示 int i; char str220; cout str2; ofstream fout (str2); if (!fout) cout n不能打开这个文件 str2 endl; exit(1); fout endl; cout endl; for (i=0; in; i+) fout xi ; cout xi ; fout endl; cout endl; fout.close (); void main () /主函数 gsdl c(3); c.input (); /从文件读入系数矩阵A、常数向量B以及e c.a_gsdl (); c.output (); /输出结果到文件并显示 实验结果：实验三 选主元高斯消去法实验说明：消元是最简单的直接法，并且也十分有效的，列主元高斯消去法对求解一般的线性方程组都适用，同时可以用来求矩阵对应的行列式。列主元高斯消去法：列主元消元 回代实验代码： #include #include #include using namespace std; class gaus private: int n; double *a, *b; public: gaus (int nn) int i; n = nn; a = new double*n; /动态分配内存空间 for (i=0; in; i+) ai = new doublen; b = new doublen; void input (); /从文件读入系数矩阵A以及常数向量B void gauss (); void output (); /输出结果到文件并显示 gaus () int i; for (i=0; in; i+) delete ai; delete a; delete b; ; void gaus:input () /从文件读入系数矩阵A以及常数向量B int i, j; char str120; cout str1; ifstream fin (str1); if (!fin) cout n不能打开这个文件 str1 endl; exit(1); for (i=0; in; i+) /读入矩阵A for (j=0; jaij; for (i=0; ibi; /读入常数向量B fin.close (); void gaus:gauss () int *js,l,k,i,j,is; double d,t; js = new intn; l=1; for (k=0; k=n-2; k+) d=0.0; for (i=k;i=n-1;i+) for (j=k;jd) d=t; jsk=j; is=i; if (d+1.0=1.0) l=0; else if (jsk!=k) for (i=0;i=n-1;i+) t=aik; aik=aijsk; aijsk=t; if (is!=k) for (j=k;j=n-1;j+) t=akj; akj=aisj; aisj=t; t=bk; bk=bis; bis=t; if (l=0) delete js; cout n系数矩阵奇异！无解. endl; return; d=akk; for (j=k+1;j=n-1;j+) akj=akj/d; bk=bk/d; for (i=k+1;i=n-1;i+) for (j=k+1;j=n-1;j+) aij=aij-aik*akj; bi=bi-aik*bk; d=an-1n-1; if (fabs(d)+1.0=1.0) delete js; cout n系数矩阵奇异！无解. =0;i-) t=0.0; for (j=i+1;j=0;k-) if (jsk!=k) t=bk; bk=bjsk; bjsk=t; delete js; void gaus:output () /输出结果到文件并显示 int i; char str220; cout str2; ofstream fout (str2); if (!fout) cout n不能打开这个文件 str2 endl; exit(1); fout endl; cout endl; for (i=0; in; i+) fout bi ; cout bi ; fout endl; cout endl; fout.close (); void main () gaus c(3); c.input (); /从文件读系数矩阵A以及常数向量B c.gauss (); c.output (); /输出结果到文件并显示 实验结果：实验四 样条插值实验说明：实现样条插值。实验代码： #include #include #include #include using namespace std; class splin3 private: int n, m; double *x, *y, *dy, *ddy, *t, *z, *dz, *ddz, integ; double *s; public: splin3 (int nn, int mm) n = nn; m = mm; x = new doublen; /动态分配内存 y = new doublen; dy = new doublen; ddy = new doublen; s = new doublen; t = new doublem; z = new doublem; dz = new doublem; ddz = new doublem; void input (); /由文件读入n个数据点(x, y)以及m个指定插值点t void interp (); /执行三次样条函数插值、微商与积分 void output (); /输出n个数据点及其一阶与二阶导数值m个插值点t处的函数值及其一阶与二阶导数值与积分值到文件并显示 splin3 () delete x, y, dy, ddy, t, z, dz, ddz, s; ; void splin3:input () /由文件读入n个数据点(x, y)以及m个指定插值点t int k; ifstream fin (ain.txt); if (!fin) cout n不能打开文件 resread.txtendl; exit(1); for (k=0; kxk; fin yk; for (k=0; ktk; /读入m个插值点 fin.close (); void splin3:interp () /执行三次样条函数插值、微商与积分 int i,j; double h0,y0,h1,y1,alpha,beta,u; h0=xn-1-xn-2; y0=yn-1-yn-2; dy0=0.0; ddy0=0.0; ddyn-1=0.0; s0=1.0; sn-1=1.0; for (j=1;j=n-1;j+) h1=h0; y1=y0; h0=xj-xj-1; y0=yj-yj-1; alpha=h1/(h1+h0); beta=3.0*(1.0-alpha)*y1/h1+alpha*y0/h0); if (j=1;j-) sj=dyj*sj+1+sj; ddyj=dyj*ddyj+1+ddyj; dyn-2=(beta-alpha*ddy1-(1.0-alpha)*ddyn-2)/ (alpha*s1+(1.0-alpha)*sn-2+2.0); for (j=2;j=n-1;j+) dyj-2=sj-1*dyn-2+ddyj-1; dyn-1=dy0; for (j=0;j=n-2;j+) sj=xj+1-xj; for (j=0;j=n-2;j+) h1=sj*sj; ddyj=6.0*(yj+1-yj)/h1-2.0*(2.0*dyj+dyj+1)/sj; h1=sn-2*sn-2; ddyn-1=6.*(yn-2-yn-1)/h1+2.*(2.*dyn-1+dyn-2)/sn-2; integ=0.0; for (i=0;i=n-2;i+) h1=0.5*si*(yi+yi+1); h1=h1-si*si*si*(ddyi+ddyi+1)/24.0; for (j=0;j=xn-1) h0=h0-(xn-1-x0); while (h0xi+1) i=i+1; u=h0; h1=(xi+1-u)/si; h0=h1*h1; zj=(3.0*h0-2.0*h0*h1)*yi; zj=zj+si*(h0-h0*h1)*dyi; dzj=6.0*(h0-h1)*yi/si; dzj=dzj+(3.0*h0-2.0*h1)*dyi; ddzj=(6.0-12.0*h1)*yi/(si*si); ddzj=ddzj+(2.0-6.0*h1)*dyi/si; h1=(u-xi)/si; h0=h1*h1; zj=zj+(3.0*h0-2.0*h0*h1)*yi+1; zj=zj-si*(h0-h0*h1)*dyi+1; dzj=dzj-6.0*(h0-h1)*yi+1/si; dzj=dzj+(3.0*h0-2.0*h1)*dyi+1; ddzj=ddzj+(6.0-12.0*h1)*yi+1/(si*si); ddzj=ddzj-(2.0-6.0*h1)*dyi+1/si; void splin3:output () /输出n个数据点及其一阶与二阶导数值m个插值点t处的函数值及其一阶与二阶导数值与积分值到文件并显示 int i; ofstream fout (aout.txt); if (!fout) cout n不能打开文件 aout.txtendl; exit(1); cout setw(4)xi setw(6)yiendl; for (i=0; in; i+) fout setw(3)xi setw(12)yiendl; cout setw(3)xi setw(12)yiendl; fout endl; cout endl; cout setw(4)ti setw(6)ziendl; for (i=0; im; i+) fout setw(3)ti setw(12)ziendl; cout setw(3)ti setw(12)ziendl; fout endl; cout endl; fout.close (); void main () /主函数 int i; splin3 solution(6,10); solution.input (); /由文件读入n个数据点(x, y)以及m个指定插值点t erp (); /执行三次样条函数插值、微商与积分 solution.output (); /输出n个数据点及其一阶与二阶导数值m个插值点t处的函数值及其一阶与二阶导数值与积分值到文件并显示 实验结果：实验五 最小二乘法实验说明：除了插值，逼近复杂函数f(x)的另一种方法是拟合。结点数据比较多，或者不够精确，往往采用拟合的方法。用来拟合的简单函数p(x)，不需严格通过给定结点，但在这些结点处总体误差达到极小，即满足最小二乘拟合的条件。线性拟合是最简单的拟合方法，许多非线性问题都可以转化为线性拟合。令线性函数为a0+a1x，要确定a0和a1，最终求解如下线性方程组：实验代码：#include #include #include using namespace std;class pir1private: int n;double *x, *y,a,b;public:pir1 (int nn)n = nn;x = new doublen; /动态分配内存y = new doublen;void input (); /由文件读入n个数据点(x, y)void fit (); /执行最小二乘直线拟合void output (); /输出拟合直线pir1 () delete x, y; ; void pir1:input () int k;ifstream fin (bin.txt);if (!fin) cout n不能打开文件 bin.txtendl; exit(1); for (k=0; kxk;fin yk; fin.close (); void pir1:fit () int k;double s1,s2,s3,s4;for(s1=0,s2=0,s3=0,s4=0,k=0;kn;k+)s1=s1+xk;s2=s2+yk;s3=s3+(xk*xk);s4=s4+(xk*yk);b=(s4/s1*double(n)-s2)/(s3/s1*double(n)-s1);a=(s2-(s1*b)/double(n);void pir1:output ()int i;ofstream fout (bout.txt);if (!fout) cout n不能打开文件 bout.txtendl; exit(1); fout a b0)cout 拟合曲线为y=a+bxendl;elseif(b=0)cout 拟合曲线为y=aendl;elsecout 拟合曲线为y=abxendl;fout.close ();void main () /主函数pir1 solution(5);solution.input();solution.fit();solution.output();实验结果：实验六 龙贝格算法实验说明：在许多实际问题中，常常需要计算定积分的值。根据微积分学基本定理，若被积函数f(x)在区间a,b上连续，只要能找到f(x)的一个原函数F(x)，便可利用牛顿-莱布尼兹公式求得积分值。但是在实际使用中，往往遇到如下困难，而不能使用牛顿-莱布尼兹公式。（1） 找不到用初等函数表示的原函数（2） 虽然找到了原函数，但因表达式过于复杂而不便计算（3） f(x)是由测量或计算得到的表格函数由于以上种种困难，有必要研究积分的数值计算问题。利用插值多项式 则积分转化为，显然易算。称为插值型求积公式。最简单的插值型求积公式是梯形公式和Simpson公式，。当求积结点提供较多，可以分段使用少结点的梯形公式和Simpson公式，并称为复化梯形公式、复化Simpson公式。如步长未知，可以通过误差限的控制用区间逐次分半的策略自动选取步长的方法称自适应算法。梯形递推公式给出了区间分半前后的递推关系。由梯形递推公式求得梯形序列，相邻序列值作线性组合得Simpson序列, Simpson序列作线性组合得柯特斯序列, 柯特斯序列作线性组合的龙贝格序列。若|R2-R1|e,则输出R2;否则依此类推。如此加工数据的过程叫龙贝格算法，如下图所示：图2.1龙贝格算法数据加工流程龙贝格算法大大加快了误差收敛的速度，由梯形序列O(h2) 提高到龙贝格序列的O(h8)实验代码： #include #include #include using namespace std; class romb private: double a, b, eps, integ; public: romb (double aa, double bb, double es)/顺序提供a,b,eps值的构造函数 a = aa; b = bb; eps = es; void integration (); void output (); double func (double); ; void romb:integration () /执行Romberg求积法 int m,n,i,k; double y10,h,ep,p,x,s,q; h=b-a; y0=h*(func (a)+func (b)/2.0; m=1; n=1; ep=eps+1.0; while (ep=eps)&(m=9) p=0.0; for (i=0;i=n-1;i+) x=a+(i+0.5)*h; p=p+func (x); p=(y0+h*p)/2.0; s=1.0; for (k=1;k=m;k+) s=4.0*s; q=(s*p-yk-1)/(s-1.0); yk-1=p; p=q; ep=fabs(q-ym-1); m=m+1; ym-1=q; n=n+n; h=h/2.0; integ = q; void romb:output () /输出积分值到文件并显示 char str220; cout str2; ofstream fout (str2); if (!fout) cout n不能打开这个文件 str2 endl; exit(1); fout integ endl; cout 积分值 = integ endl; fout.close (); double romb:func (double x) if (x=0) return 1; elsereturn (double)(sin(x)/x); void main () /主函数 romb solution(0.0, 1.0, 0.000001); /创建对象并顺序提供a, b, eps值 egration (); solution.output (); 实验结果：实验七 矩阵三角分解的LU算法实验说明：算法原理：应用高斯消去法解阶线性方程经过步消去后得出一个等价的上三角形方程组，对上三角形方程组用逐步回代就可以求出解来。这个过程也可通过矩阵分解来实现。将非奇异阵分解成一个下三角阵和上三角阵的乘积称为对矩阵的三角分解，又称分解。根据分解,将分解为形式,简化了求解问题。程序框图：变量说明：为系数矩阵元素，为常数矩阵系数，分别为下、上三角矩阵元素。实验代码： #include #include #include using namespace std; class lluu private: int n; double *a, *l, *u; public: lluu (int nn) int i; n = nn; a = new double*n; /动态分配内存空间 for (i=0; in; i+) ai = new doublen; l = new double*n; for (i=0; in; i+) li = new doublen; u = new double*n; for (i=0; in; i+) ui = new doublen; void input (); /由文件读入矩阵A的元素 void lu (); /LU分解 void output (); /输出L矩阵与U矩阵到文件并显示 lluu () int i; for (i=0; in; i+) delete ai; delete a; for (i=0; in; i+) delete li; delete l; for (i=0; in; i+) delete ui; delete u; ; void lluu:input () /由文件读入矩阵A的元素 int i, j; char str120; cout str1; ifstream fin (str1); if (!fin) cout n不能打开这个文件 str1 endl; exit(1); for (i