典型数值算法的 C++语言程序设计.doc

数值计算课程设计数值方法课程设计说明书题目： 典型数值算法的C+语言程序设计 学生姓名： 晏瑞 学 号： 200912010127 院 （系）： 理学院 专 业： 数学与应用数学091班指导教师： 刘海峰 2011 年 6 月 15 日- 43 -陕 西 科 技 大 学数值计算课程设计任务书理学院 应用数学专业 数学091班级 学生： 晏瑞 题目：典型数值算法的C+语言程序设计 课程设计从 2011 年 5 月 20日起到 2011 年 6 月 25 日1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）：每人需作10个算法的程序、必做6题、自选4题。对每个算法要求用C+语言进行编程。必选题：1、经典四阶龙格库塔法解一阶微分方程组2、高斯列主元法解线性方程组3、牛顿法解非线性方程组4、龙贝格求积分算法5、三次样条插值算法（压紧样条）用C+语言进行编程计算 依据计算结果，用Matlab画图并观察三次样条插值效果。6、M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果。自选题：自选4道其他数值算法题目.每道题目重选次数不得超过5次.2、对课程设计成果的要求包括图表、实物等硬件要求：1）提交课程设计报告按照算法要求，用C+语言设计和开发应用程序，提交由算法说明；程序设计说明；系统技术文档 (包括系统各模块主要流程图，软件测试方案与测试记录、软件调试和修改记录、测试结论、运行情况记录)，系统使用说明书，源程序代码为附录构成的课程设计报告。2）课程设计报告版式要求打印版面要求：A4纸，页边距：上2cm，下2cm，左2.5cm、右2cm；字体：正文宋体、小四号；行距：固定值20；页眉1.5cm ，页脚1.75cm；页码位于页脚居中打印；奇数页页眉“数值计算课程设计”，偶数页页眉“算法名称”，页眉宋体小5号；段落及层次要求：每节标题以四号黑体左起打印（段前段后各0.5行），节下为小节，以小四号黑体左起打印（段前段后各0.5行）。换行后以小四号宋体打印正文。节、小节分别以1、1.1、1.1.1依次标出，空一字符后接各部分的标题。当论文结构复杂，小节以下的标题，左起顶格书写，编号依次用（1）、（2）或1）、2）顺序表示。字体为小四号宋体。 对条文内容采用分行并叙时，其编号用（a）、（b）或a）、b）顺序表示，如果编号及其后内容新起一个段落，则编号前空两个中文字符。3）设计报告装订顺序与规范封面数值计算课程设计任务书目录数值计算设计课程设计报告正文设计体会及今后的改进意见参考文献（资料）左边缘装订3、课程设计工作进度计划：时间设计任务及要求第16周编写和调试程序并按要求撰写设计报告 指导教师： 日期： 教研室主任： 日期： 目录1.经典四阶龙格库塔法解一阶微分方程- 5 -1.1算法说明：- 5 - 1.2，算法程序：- 6 -1.3，运行结果- 8 -2.高斯列主元法解线性方程组- 9 -2.1 算法说明：- 9 -2.2 算法程序：- 9 -2.3 运行结果：- 13 -3.牛顿法解非线性方程组- 14 -3.1算法说明- 14 -3.2 算法程序：- 14 -3.3运行结果：- 19 -4.龙贝格求积分算法- 20 -4.1算法说明- 20 -4.2 算法程序：- 20 -4.3 运行结果：- 22 -5.三次紧压样条插值- 23 -5.1 算法说明：- 23 -5.2 算法程序：- 23 -5.3 运行结果：- 25 -6.M次多项式曲线拟合- 26 -6.1算法说明- 26 -6.2 算法程序：- 26 -7.不动点法解非线性方程- 33 -7.1算法说明- 33 -7.3 运行结果：- 35 -8.二分法解非线性方程- 36 -8.1 算法说明：- 36 -8.4运行结果：- 38 -9，龙格-库塔法解微分方程- 39 -9.1 算法说明：- 39 -9.2 算法程序：- 39 -9.3运行界面：- 41 -10，递归梯形公式- 42 -10.1，算法说明：- 42 -10.2，算法程序：- 42 -10.3 运行结果：- 43 -1.经典四阶龙格库塔法解一阶微分方程1.1算法说明： 龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。对于一精度的欧拉公式有： yi+1=yi+h*K1 K1=f(xi,yi) 当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进欧拉公式： yi+1=yi+h*( K1+ K2)/2 K1=f(xi,yi) K2=f(xi+h,yi+h*K1) 依次类推，如果在区间xi,xi+1内多预估几个点上的斜率值K1、K2、Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格库塔公式，也就是在工程中应用广泛的经典龙格库塔算法： yi+1=yi+h*( K1+ 2*K2 +2*K3+ K4)/6 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) 通常所说的龙格-库塔法是指四阶而言的，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式 1.2，算法程序：#include#include#define M 10using namespace std;int main()long double feval(long double ,long double ,long double );long double feval1(long double ,long double ,long double );long double f4,g4;long double h,a,b,xa,ya;long double xM+1,yM+1,TM+1;coutab;coutendl;/xa,ya是初值coutxaya;coutendl;/h是步长h=(b-a)/M;x0=xa;y0=ya;for(int i=0;i=M;i+)Ti=a+h*i; /给T赋值for(int k=0;k=M;k+)f0=feval(Tk,xk,yk); g0=feval1(Tk,xk,yk);f1=feval(Tk+h/2,xk+h/2*f0,yk+h/2*g0);g1=feval1(Tk+h/2,xk+h/2*f0,yk+h/2*g0);f2=feval(Tk+h/2,xk+h/2*f1,yk+h/2*g1);g2=feval1(Tk+h/2,xk+h/2*f1,yk+h/2*g1);f3=feval(Tk+h,xk+h*f2,yk+h*g2);g3=feval1(Tk+h,xk+h*f2,yk+h*g2);xk+1=xk+h/6*(f0+2*f1+2*f2+f3);yk+1=yk+h/6*(g0+2*g1+2*g2+f3);coutsetw(12)Tsetw(12)xsetw(12)yendl;for(k=0;k=M;k+)coutsetw(12)Tksetw(12)xksetw(12)ykendl;return 0;long double feval(long double t,long double x,long double y)long double f;f=x+2*y;return f; long double feval1(long double t,long double x,long double y)long double f;f=3*x+2*y;return f;1.3，运行结果 2.高斯列主元法解线性方程组 2.1 算法说明：首先将线性方程组做成增广矩阵，对增广矩阵进行行变换。对于元素，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该最大元素所在的行与第i行交换，然后采用高斯消元法用新得到的所在的第i 行消去第i行以下的元素。依次进行直到。从而得到上三角矩阵。2.2 算法程序：本程序包含enter.h，remove.h，judge.h及main.cpp这四个函数1，main.cpp函数#include#include#include enter.h /包含名为enter的头函数#include remove.h /包含名为remove的头函数 #include judge.h /包含名为judge的头函数int main()int row,col;coutrow;coutcol;double a100100,*p1,*p100;for(int i=0;irow;i+)pi=ai;p1=p;enter(p1,row,col); /调用enter函数remove(p1,row,col); /调用remove函数int number=judge(p1,row,col); /调用judge函数 if(number=0)cout-结果-endl;cout此方程组无解endl; /无解时予以提示 if(number=-1)cout-结果-endl;cout此方程组有无穷多解endl; /有无穷多解时予以提示if(number=1)cout-结果-endl;cout此方程组有唯一解=0;i-)double t=aicol-1;for(int j=col-2;ji;j-)t=t-aij*xj;xi=t/aii;for(i=0;icol-1;i+)coutxi+1=xit; /依次输出方程组的解coutendl;return 0;2，enter.h函数 /输入系数矩阵void enter(double *q,int m,int n)cout请按行输入未知数前面的系数和等式右边的常数：endl;for(int i=0;im;i+)for(int j=0;jqij;3，remove.h函数/寻找列主元，并移动该行，最后化为上三角矩阵void remove(double *q,int m,int n) int min=m;if(n-1m)min=n-1;for(int i=0;imin-1;i+)int k=i;double max=qii;for(int j=i+1;jfabs(max)max=qji;k=j; /找到第i列从aii开始的绝对值最大的元素if(k!=i)for(int j=0;jn;j+)double mat=qij;qij=qkj;qkj=mat; /通过换行以保证主对角线上的元素是其所在位置及以下元素中绝对值最大的一个int t=0;for(j=i;jm;j+)if(qji=0)t+; /在消元前判断aii及其所在列以下元素是否都为零，不都为零再进行消元if(t!=m-i) for(int j=i+1;j=i;k-) qjk=qjk-qik*qji/qii;for(i=0;in-2;i+)for(int j=i+1;jm;j+)qji=0; /保证经消元后得到上三角阵4，judge.h函数/用增广矩阵的秩判断该矩阵所对应的方程组的解的情况int judge(double *q,int m,int n) int r1=m;for(int i=0;im;i+)int g=0;for(int j=0;jn;j+)if(qij=0)g+;if(g=n)r1=r1-1; /求出增广矩阵的秩int r2=m;for(i=0;im;i+)int h=0;for(int j=0;jn-1;j+)if(qij=0)h+;if(h=n-1)r2=r2-1; /求出系数矩阵的秩if(r1=n-1&r2=n-1)return 1;else if(r1=r2&r1n-1)return -1; /据不同的秩的情况返回不同的值 else return 0;2.3 运行结果：3.牛顿法解非线性方程组3.1算法说明设已知。第1步：计算函数第2步：计算雅可比矩阵第3步：求线性方程组的解。第4步：计算下一点重复上述过程。3.2 算法程序：#include#include#define N 3 / 非线性方程组中方程个数、未知量个数 #define Epsilon 0.0001 / 差向量1范数的上限#define Max 3 /最大迭代次数using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN);/计算向量函数的因变量向量yyNvoid ffjacobian(float xxN,float yyNN);/计算雅克比矩阵yyNNvoid inv_jacobian(float yyNN,float invNN);/计算雅克比矩阵的逆矩阵invvoid newdundiedai(float x0N, float invNN,float y0N,float x1N);/由近似解向量 x0 计算近似解向量 x1float x0N=0,0,0,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,j,iter=0;/如果取消对x0的初始化，撤销下面两行的注释符，就可以由键盘向x0读入初始近似解向量/for( i=0;ix0i;cout初始近似解向量：endl;for (i=0;iN;i+) coutx0i ; coutendl;coutendl;doiter=iter+1;cout第 iter 次迭代开始endl;/计算向量函数的因变量向量 y0ff(x0,y0);/计算雅克比矩阵 jacobianffjacobian(x0,jacobian);/计算雅克比矩阵的逆矩阵 invjacobianinv_jacobian(jacobian,invjacobian);/由近似解向量 x0 计算近似解向量 x1 newdundiedai(x0, invjacobian,y0,x1);/计算差向量的1范数errornormerrornorm=0;for (i=0;iN;i+)errornorm=errornorm+fabs(x1i-x0i);if (errornormEpsilon) break;for (i=0;iN;i+)x0i=x1i; while (iterMax);return 0;void ff(float xxN,float yyN)float x,y,z; int i;x=xx0;y=xx1;z=xx2; yy0=x*x-x+y*y+z*z-5;yy1=x*x+y*y-y+z*z-4;yy2=x*x+y*y+z*z+z-6; cout向量函数的因变量向量是： endl;for( i=0;iN;i+)coutyyi ; coutendl;coutendl;void ffjacobian(float xxN,float yyNN) float x,y,z; int i,j; x=xx0;y=xx1; z=xx2;/jacobian have n*n elementyy00=2*x-1;yy01=2*y;yy02=2*z;yy10=2*x;yy11=2*y-1;yy12=2*z;yy20=2*x;yy21=2*y;yy22=2*z+1;cout雅克比矩阵是： endl;for( i=0;iN;i+) for(j=0;jN;j+) coutyyij ; coutendl; coutendl;void inv_jacobian(float yyNN,float invNN)float augNN2,L; int i,j,k;cout开始计算雅克比矩阵的逆矩阵 ：endl; for (i=0;iN;i+) for(j=0;jN;j+) augij=yyij; for(j=N;jN2;j+)if(j=i+N) augij=1;else augij=0;for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl;coutendl;for (i=0;iN;i+) for (k=i+1;kN;k+) L=-augki/augii;for(j=i;jN2;j+) augkj=augkj+L*augij; for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl;cout0;i-) for (k=i-1;k=0;k-) L=-augki/augii;for(j=N2-1;j=0;j-) augkj=augkj+L*augij; for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl;cout=0;i-) for(j=N2-1;j=0;j-)augij=augij/augii;for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij ; coutendl; for(j=N;jN2;j+) invij-N=augij;coutendl;cout雅克比矩阵的逆矩阵： endl;for (i=0;iN;i+) for(j=0;jN;j+) coutinvij ; coutendl;coutendl;void newdundiedai(float x0N, float invNN,float y0N,float x1N)int i,j;float sum=0;for(i=0;iN;i+) sum=0; for(j=0;jN;j+) sum=sum+invij*y0j;x1i=x0i-sum;cout近似解向量：endl;for (i=0;iN;i+) coutx1i ;coutendl;coutendl;3.3运行结果：结果为：x=-1.36628;y=-0.366281;z=2.366284.龙贝格求积分算法4.1算法说明生成的逼近表，并以为最终解来逼近积分逼近存在于一个特别的下三角矩阵中，第0列元素用基于个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为当时，程序在第行结束。4.2 算法程序：/龙贝格方法求定积分/函数的修改可以在feval函数的定义中修改#include#include#includeusing namespace std;#define N 4 /定义全局变量Nint main()long double feval(long double );/声明在指定函数的指定点处的函数long double a,b,h;long double R1010=0;long double x;int M=1;long double err=1;int J=0;long double s=0;coutab;h=b-a;R00=h*(feval(a)+feval(b)/2;while(err0.0001)&(JN)|(J4) /函数的循环J=J+1;h=h/2;s=0;for(int p=1;p=M;p+)x=a+h*(2*p-1); s=s+feval(x);RJ0=RJ-10/2+h*s;M=2*M;for(int K=1;K=J;K+)RJK=RJK-1+(RJK-1-RJ-1K-1)/(4K-1);err=fabs(RJ-1J-1-RJK);cout该矩阵为：endl;for(int i=0;i=J;i+)for(int j=0;j=J;j+)coutsetiosflags(ios:fixed)setprecision(8)setw(12)Rij;coutendl;cout该式子在a,b内的积分为RJJ;coutendl;return 0;long double feval(long double x) /函数的定义long double f;f=(x*x+x+1)*cos(x);return(f);4.3 运行结果：5.三次紧压样条插值5.1 算法说明：三次紧压样条，确定，（如果导数已知，这是“最佳选择”），。5.2 算法程序：#include#includeusing namespace std;int main()const int n=4; const int N=n-1;/3long double Xn=0,1,2,3,Yn=0,0.5,2.0,1.5,dx0=0.2,dxn=-1,Hn-1,H1n-1,Dn-1,Un-2,AN-2,BN-1,CN-1,MN+1,temp,SN4;cout求解经过点(0，0.0),(1，0.5),(2，2.0)和(3，1.5),而且一阶导数边界条件S(0)=0.2和S(3)=-1的三次压紧样条曲线:;for(int i=0;in-1;i+)Hi=Xi+1-Xi;for(int i1=0;i1n-1;i1+)H1i1=Yi1+1-Yi1;for(int i2=0;i2n-1;i2+)Di2=H1i2/Hi2;for(int i3=0;i3n-1;i3+)Ui3=6*(Di3+1-Di3);for(int i4=0;i4N-2;i4+)Ai4=Hi4+1;for(int i5=0;i5N-1;i5+)Bi5=2*(Hi5+Hi5+1);for(int i6=0;i6N-1;i6+)Ci6=Hi6+1;B0=B0-H0/2;U0=U0-3*(D0-dx0);BN-2=BN-2-HN-1/2;UN-2=UN-2-3*(dxn-DN-1);for(int k=1;k0;k1-)Mk1=(Uk1-1-Ck1-1*Mk1+1)/Bk1-1;M0=3*(D0-dx0)/H0-M1/2;MN=3*(dxn-DN-1)/HN-1-MN-1/2;for(int k2=0;k2N;k2+)Sk20=(Mk2+1-Mk2)/(6*Hk2);Sk21=Mk2/2;Sk22=Dk2-Hk2*(2*Mk2+Mk2+1)/6;Sk23=Yk2;coutS=:endl;for(int j=0;jN;j+)for(int j1=0;j14;j1+)coutsetw(12)Sjj1;coutendl;return 0;5.3 运行结果：6.M次多项式曲线拟合6.1算法说明设有N个点，横坐标是确定的。最小二乘抛物线的系数表示为求解A，B和C的线性方程组为6.2 算法程序：#include#include#includeusing namespace std;const int n=6,m=5,N=6,N2=2*N;int main()long double inv_jacobian(long double yyNN,long double invNN);long double Xn=0.25,1.0,1.5,2.0,2.4,5.0,Yn=23.1,1.68,1.0,0.84,0.826,1.2576,Bm+1,Fnm+1=0,F1m+1n,F2nm+1,Am+1m+1,Y1n,invm+1m+1,Cm+1,B1m+1;for(int i=0;im+1;i+)for(int j=0;jn;j+)Fji=pow(Xj,i);/F的转置for(int i1=0;i1m+1;i1+)for(int j1=0;j1n;j1+)F1i1j1=Fj1i1;/输出F1coutF的转置F=:endl;for(int p1=0;p1m+1;p1+)for(int k1=0;k1n;k1+)coutsetw(13)F1p1k1;coutendl;coutendl;/输出FcoutF=:endl;for(int p2=0;p2n;p2+)for(int k2=0;k2m+1;k2+)coutsetw(13)Fp2k2;coutendl;coutendl;coutendl;/F的转置F1乘以F得A矩阵for(int i2=0;i2m+1;i2+) for(int j2=0;j2n;j2+) for(int k1=0;k1m+1;k1+) F2j2k1=F1i2j2*Fj2k1;/每次输出F2cout第i2+1次乘法：endl;for(int p3=0;p3n;p3+)for(int k3=0;k3m+1;k3+)coutsetw(13)F2p3k3;coutendl;coutendl;/得到矩阵A for(int j3=0;j3m+1;j3+) Ai2j3=0;/必须给此数赋初值 for(int i3=0;i3n;i3+) Ai2j3=Ai2j3+F2i3j3;/输出AcoutA=F*F=:endl;for(int p=0;pm+1;p+)for(int k=0;km+1;k+)coutsetw(13)Apk;coutendl;coutendl;for(int i4=0;i4m+1;i4+) for(int j4=0;j4n;j4+) Y1j4=F1i4j4*Yj4;/每次输出Y1cout第i4+1次乘法：endl; for(int p4=0;p4n;p4+) coutsetw(13)Y1p4;coutendl;/得到矩阵BBi4=0;/必须给此数赋初值 for(int i5=0;i5m+1;i5+) Bi4=Bi4+Y1i5;coutB=F*Y=:endl;for(int i6=0;i6n;i6+)coutsetw(13)Bi6;coutendl;invm+1m+1=inv_jacobian(A,inv);coutA矩阵的逆矩阵： endl;for (i=0;iN;i+) for(int j0=0;j0N;j0+) coutsetw(13)invij0; coutendl;coutendl;for(int qq=0;qqm+1;qq+) for(int w=0;wn;w+) B1w=invqqw*Bw;/每次输出B1cout第qq+1次乘法：endl; for(int ww=0;wwn;ww+) coutsetw(13)B1ww;coutendl;/得到矩阵CCqq=0;/必须给此数赋初值 for(int e=0;em+1;e+) Cqq=Cqq+B1e;coutC=A*B=endl;for(int tt=0;ttm+1;tt+)coutsetw(13)Ctt;coutendl;return 0;long double inv_jacobian(long double yyNN,long double invNN)long double augNN2,L; int i,j,k; cout开始计算A矩阵的逆矩阵 ：endl; for (i=0;iN;i+) for(j=0;jN;j+) augij=yyij; for(j=N;jN2;j+)if(j=i+N) augij=1;else augij=0; cout第零次消去：endl; for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij/; coutendl;coutendl;for (i=0;iN;i+) for (k=i+1;kN;k+) L=-augki/augii;for(j=i;jN2;j+) augkj=augkj+L*augij; cout第一次消去：endl; for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij/; coutendl;cout0;i-) for (k=i-1;k=0;k-) L=-augki/augii;for(j=N2-1;j=0;j-) augkj=augkj+L*augij; cout第二次消去：endl;for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij/; coutendl;cout=0;i-) for(j=N2-1;j=0;j-)augij=augij/augii;cout第三次消去：endl;for (i=0;iN;i+) for(j=0;jN2;j+) coutaugij/; coutendl; for(j=N;jN2;j+) invij-N=augij;coutendl;return(invm+1m+1);6.3 运行结果：7.不动点法解非线性方程7.1算法说明先将改写成然后对进行迭代，即 其中然后判断是否成立，成立则返回，不成立就重复以上步骤7.2 算法程序：/不动点法求非线性方程#include#include#define N 10using namespace std;int main()double