典型数值算法的说明 C++语言程序的设计开发.docx

1、典型数值算法的说明 C+语言程序的设计开发数值方法课程设计说明书题目： 典型数值算法的C+语言程序设计 *学 号： * 院 （系）： 理学院 专 业： 数学与应用数学091班*2011 年 6 月 15 日陕 西 科 技 大 学数值计算课程设计任务书理学院 应用数学专业 数学091班级 学生： 晏瑞 题目：典型数值算法的C+语言程序设计 课程设计从 2011 年 5 月 20日起到 2011 年 6 月 25 日1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）：每人需作10个算法的程序、必做6题、自选4题。对每个算法要求用C+语言进行编程。必选题：1、经典四阶龙格库塔法解一阶微分

2、方程组2、高斯列主元法解线性方程组3、牛顿法解非线性方程组4、龙贝格求积分算法5、三次样条插值算法（压紧样条）用C+语言进行编程计算 依据计算结果，用Matlab画图并观察三次样条插值效果。6、M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果。自选题：自选4道其他数值算法题目.每道题目重选次数不得超过5次.2、对课程设计成果的要求包括图表、实物等硬件要求：1）提交课程设计报告按照算法要求，用C+语言设计和开发应用程序，提交由算法说明；程序设计说明；系统技术文档 (包括系统各模块主要流程图，软件测试方案与测试记录、软件调试和修改记录、测试结论、运行情况记录)，系统使用说明书，源程

3、序代码为附录构成的课程设计报告。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）

4、、2）顺序表示。字体为小四号宋体。 对条文内容采用分行并叙时，其编号用（a）、（b）或a）、b）顺序表示，如果编号及其后内容新起一个段落，则编号前空两个中文字符。3）设计报告装订顺序与规范封面数值计算课程设计任务书数值计算设计课程设计报告正文设计体会及今后的改进意见参考文献（资料）左边缘装订3、课程设计工作进度计划：时间设计任务及要求第16周编写和调试程序并按要求撰写设计报告 指导教师： 日期： 教研室主任： 日期： 1.经典四阶龙格库塔法解一阶微分方程1.1算法说明： 龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所

5、以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。对于一精度的欧拉公式有： 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*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格库塔公式，也就是在工程中应用广泛

6、的经典龙格库塔算法： 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 );l

7、ong 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

8、(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(1

9、2)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 算法说明：首先将线性方程组做成

10、增广矩阵，对增广矩阵进行行变换。对于元素，在第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 ma

11、in() 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; /有

12、无穷多解时予以提示 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;

13、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; /通过换行以保证主对角线上的元素是其所在位置及以

14、下元素中绝对值最大的一个 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+) in

15、t 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步：计算函数 第

16、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);/计算雅克比矩阵yyNNvoi

17、d 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+

18、) coutx0i ; coutendl;coutendl;do iter=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+) e

19、rrornorm=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 ; couten

20、dl; coutendl;void ffjacobian(float xxN,float yyNN) float x,y,z; int i,j; x=xx0; y=xx1; z=xx2; /jacobian have n*n element yy00=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

21、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;

22、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+

23、) 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算法说明 生成的逼近表，并以为最终解来逼近积分