数值计算-课程设计选编

1、数值计算课程设计说明书题目： 典型数值算法的C+语言程序设计 院 系： 理学院 专业班级： 信息131 学 号： 202112030120 学生姓名： 郝腾宇 指导教师： 刘海峰 2021年 7月 10 日陕 西 科 技 大 学数值计算课程设计任务书理学院 信息与计算科学 专业 信息131 班级 学生： 郝腾宇 题目：典型数值算法的C+语言程序设计 课程设计从 2021 年 5 月 20 日起到 2021 年 7月 10 日1、课程设计的内容和要求包括原始数据、技术要求、工作要求等：每人需作10个算法的程序、必做6题、自选4题。对每个算法要求用C+语言进行编程。必选题：1、高斯列主元法解线性方

2、程组2、牛顿法解非线性方程组3、经典四阶龙格库塔法解一阶微分方程组4、三次样条插值算法压紧样条用C+语言进行编程计算 依据计算结果，用Matlab画图并观察三次样条插值效果。5、龙贝格求积分算法6、M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果。自选题：自选4道其他数值算法题目.每道题目重选次数不得超过5次.2、对课程设计成果的要求包括图表、实物等硬件要求：2.1 提交课程设计报告按照算法要求，应用C+语言设计和开发算法程序，提交由：每个算法的原理与公式说明；与算法相应的程序设计说明程序中的主要变量语义说明，变量的数据类型，数据在内存中组织和存储结构说明，各函数模块的主要流程

3、图，函数功能说明，函数的形参说明，函数的调用方法说明；程序调试与实例运行记录 (包括程序调试和修改记录、测试结论、运行结果记录)，每个算法的源程序代码编入附录构成的课程设计报告。2.2 课程设计报告版式要求目录的要求：居中打印目录二字，四号黑体，段后1行，字间空一字符；章、节、小节及其开始页码字体均为小四号宋体。节向右缩进两个字符汉字，小节及以后标题均向右缩进四个字符汉字。目录中应包含正文局部每个算法章节标题、设计体总结、无序号的“参考文献资料，目录的最后一项为哪一项“附录正文的要求：算法说明论述清楚，公式符号撰写标准，流程图图符标准， 计算正确，文字简练通顺，插图简洁标准，书写整洁。文中图、

4、表按制图要求绘制，程序调试和运行情况记录详实。打印版面要求：A4纸，页边距：上2cm，下2cmcm、右2cm；字体：正文宋体、小四号；行距：固定值20；页眉1.5cm ，页脚1.75cm；页码位于页脚居中打印；奇数页页眉“数值计算课程设计，偶数页页眉“具体算法名，页眉宋体小5号；段落及层次要求：每节标题以四号黑体左起打印段前段后各0.5行，节下为小节，以小四号黑体左起打印段前段后各0.5行。换行后以小四号宋体打印正文。章、节、小节编号分别以1、1.1、1.1.1格式依次标出，空一字符后接各局部的标题。每一章的标题都应出现在本章首页的第一行上。当课程设计报告结构复杂，小节以下的标题，左起顶格书写

5、，编号依次用1、2或1、2顺序表示。字体为小四号宋体。对条文内容采用分行并叙时，其编号用a、b或a、b顺序表示，如果编号及其后内容新起一个段落，那么编号前空两个中文字符。曲线图表要求：所有曲线、图表、线路图、流程图、程序框图、示意图等不准徒手画，必须按国家规定标准或工程要求绘制应采用计算机绘图。课程设计说明书报告中图表、公式要求如下：a图：图的名称采用中文，中文字体为五号宋体，图号图名在图片下面。引用图应在图题右上角标出文献来源。图号以章为单位顺序编号。格式为：图1-1，空一字符后接图名，比方第1章第5个图是关于高斯列主元法解方程组算法图，图的下方的图号图名应为：图1-5 高斯列主元法解方程组

6、算法图。b表格：表的名称及表内文字采用中文，中文字体为五号宋体，表号表名在表格上面。表号以章为单位顺序编号，表内必须按规定的符号标注单位。格式为：表1-1，空一字符后接表格名称。比方第4章第1个表是关于三次样条插值的插值点列表，表的上方表号表名那么应为：表4-1 插值点的列表。c公式：公式书写应在文中另起一行，居中排列。公式序号按章顺序编号。字体为五号宋体，序号靠页面右侧。比方第3章第1个公式其编号那么应为：3-1。设计体会及今后的改良意见：设计总结要写出算法理解，编程经验等技术性、学术性总结；体会要简洁、真实、深刻，切忌空话、大话，客套话和矫揉造作之词。改良意见要合理、中肯。参考文献的要求：

7、另起一页，居中打印参考文献四字四号黑体，段前段后1行，字间空一字符；另起一行，按论文中参考文献出现的先后顺序用阿拉伯数字连续编号参考文献编号应在正文中标注出；参考文献中每条工程应齐全字体均为小四号宋体。格式：编号作者.论文或著作名称.期刊名或出版社.出版时间。期刊应注明第几期、起止页数包括论著。参考文献中条目要符合科技文献引用文献条目书写的国家标准标准。2.3 设计报告装订顺序与标准封面数值计算课程设计任务书目录数值计算课程设计报告正文设计体会及今后的改良意见参考文献(无需加目录序号)附录无需加目录序号)左边缘装订3、课程设计工作进度方案：时间设计任务及要求第18周编写和调试程序并按要求撰写设

8、计报告 目 录1 高斯列主元法解线性方程组11.1 高斯列主元法解线性方程组算法思想11.2 高斯列主元法解线性方程组流程图11.3 高斯列主元法解线性方程组运行结果22 牛顿法解非线性方程组32.1 牛顿法解非线性方程组算法思想32.2 牛顿法解非线性方程组流程图42.3 牛顿法解非线性方程组运行结果52.4 牛顿法解非线性方程组手工计算63 经典四阶龙格库塔法解一阶微分方程组73.1 经典四阶龙格库塔法解一阶微分方程算法思想73.2 经典四阶龙格库塔法解一阶微分方程流程图83.3 经典四阶龙格库塔法解一阶微分方程运行结果94 三次样条插值算法压紧样条104.1 三次样条插值算法说明104.

9、2 三次样条插值算法流程图 114.3 三次样条插值运行结果124.4 Matlab画图并观察三次样条插值效果135 龙贝格求积分算法14龙贝格算法说明14龙贝格算法流程图15龙贝格算法运行结果166 M次多项式拟合176.1 M次多项式拟合算法思想176.2 M次多项式拟合流程图176.3 程序运行结果186.4 MATLAB拟合结果187 雅克比迭代法197.1 雅克比迭代法算法说明197.2 雅克比迭代法流程图197.3 雅克比迭代法程序运行结果208 高斯-赛德尔迭代法218.1 算法说明218.2 高斯-赛德尔迭代法流程图228.3 高斯-赛德尔迭代法流运行结果239 复化辛普森公式

10、求积分249.1 复化辛普森公式求积分算法说明249.2 复化辛普森公式求积分流程图259.3 复化辛普森公式运行结果2610 最小二乘线性拟合2710.1 最小二乘线性拟合算法说明2710.2 最小二乘法进行线性拟合流程图2710.3 最小二乘法进行线性拟合运行结果2811 设计体会及改良意见29参 考 文 献30附 录31附录一 高斯列主元法程序代码31附录二 牛顿法解非线性方程组的程序代码35附录三 经典四阶龙格库塔法39附录四 三次样条插值算法压紧样条41附录五 龙贝格算法程序代码44附录六 M次多项式曲线拟合程序代码46附录七 雅克比迭代程序代码50附录八 高斯赛德尔迭代程序代码52

11、附录九 复化辛普森求积公式程序代码55附录十 最小二乘法线性拟合程序代码571 高斯列主元法解线性方程组问题提出：求解线性方程组1.1 高斯列主元法解线性方程组算法思想1） 将线性方程组写成增广矩阵的形式；2） 对增广矩阵进行行变换，对元素，在第i列中及以下的元素选取绝对值最大的元素，将所有元素中最大的所在的行与第i行交换，然后采用高斯消元法使得新得到的第i行以下的元素均为零。一直重复上述过程直到。从而得到上三角矩阵。3对上三角矩阵进行回代求解，即可以得到方程组的解。1.2 高斯列主元法解线性方程组流程图输入增广矩阵Start 输入未知量的个数n高斯列主元法变换为上三角矩阵矩阵是否退化YN回代

12、求解End 图1-1 高斯列主元法流程图1.3 高斯列主元法解线性方程组运行结果图1-2 高斯列主元法程序运行界面图2 牛顿法解非线性方程组问题提出：求解非线性方程组设初值，用牛顿法计算。2.1 牛顿法解非线性方程组算法思想前提：设1计算函数2-12计算雅可比矩阵2-2(2-3)3求线性方程组的解。(2-4)4计算下一点重复上述过程，直到到达要求的精度，停止计算。2.2 牛顿法解非线性方程组流程图Start 输入精度要求和最大迭代次数数 计算计算雅可比矩阵 =- 或者精度不够Y EndN图2-1 牛顿法解非线性方程组流程图2.3 牛顿法解非线性方程组运行结果图2-2 牛顿法解非线性方程组程序运

13、行图2.4 牛顿法解非线性方程组手工计算求解非线性方程组设初值，用牛顿法计算。函数向量和雅克比矩阵为(2-5)在点2.00,0.25处的函数值为(2-6)由牛顿法可得出下面的线性方程组 解得(2-7)迭代得(2-8)手工计算与程序运行结果一致。3 经典四阶龙格库塔法解一阶微分方程组问题提出：求解线性微分方程组 3.1 经典四阶龙格库塔法解一阶微分方程算法思想龙格-库塔法的根本思想是：用在几个不同点的数值加权平均来代替的值，而使截断误差的阶数尽可能高。我们用四个不同点上的函数值的线性组合，将精度提高到四阶就可以得到四阶龙格-库塔公式。常用的标准四阶龙格-库塔公式为： (3-1)其中 (3-2)

14、经典四阶龙格库塔法解一阶微分方程流程图 Start Y N输出结果End图3-1 四阶龙格库塔法程序流程图3.3 经典四阶龙格库塔法解一阶微分方程运行结果图3-2 四阶龙格库塔法程序运行界面4 三次样条插值算法压紧样条4.1 三次样条插值算法说明设及，插值函数在每个小区间上是不超过3次的多项式且具有2阶连续导数，那么称为三阶样条插值。具体的说，三阶样条插值是满足以下条件的分段三次多项式：（1） 插值条件：(4-1)； (2)连接条件：(4-2)压紧样条公式(4-3)三阶样条插值的边界条件为1两端的一阶导数；2两端的二阶导数；3默认:自然边界条件。4.2 三次样条插值算法流程图 Start输入所

15、有点坐标及边界条件计算三次样条插值系数输出结果End图4-1 三次样条插值程序流程图4.3 三次样条插值运行结果图4-2 三次样条差值运行界面图4.4 Matlab画图并观察三次样条插值效果图4-3 三次样条插值Matlab画图5 龙贝格求积分算法问题提出:用龙贝格积分法求解以下积分：(5-1)步1 输入a,b及精度；步2 步3 对分区间a,b，并计算(5-2)步4 假设不满足终止条件，作循环：(5-3)(5-5)(5-4) 计算 对计算终止条件一般取为手工计算过程如下：开始 输入下限a、上限b、精度。 ,，计算 得到积分值结束图5-1 龙贝格积分程序流程图图5-2 龙贝格积分程序运行图6 M

16、次多项式拟合问题提出：根据4个数据点-3,30,12,14,3，求解最小二乘抛物线。6.1 M次多项式拟合算法思想6-1本算法基于最小二乘原理来通过数据点来构造M解最小二乘多项式。 以为例考虑求解最小二乘的抛物线拟合。设有个点，横坐标确实定的。最小二乘抛物线的系数表示为6-2求解的线性方程组为6-56-46-36.2 M次多项式拟合流程图解线性方程组得到拟合多项式系数Start 输入N个数据点坐标和拟合多项式系数M用x、y的对应次方构造方程 End图6-1 多项式拟合程序流程图6.3 程序运行结果图6-2 多项式拟合运行结果6.4 MATLAB拟合结果图6-3 MATLAB拟合结果7 雅克比迭

17、代法7.1 雅克比迭代法算法说明步1 取初始点，精度要求，最大迭代次数N,置步2 由计算步3 假设那么停算，输出作为方程组的近似解；步4 假设那么停算，输出迭代失败信息；否那么置，转步2.7.2 雅克比迭代法流程图开始 请输入精度要求和最大迭代次数 计算判断 否 是结束 图7-1 雅克比迭代程序流程图7.3 雅克比迭代法程序运行结果图7-2 雅克比迭代运行结果8 高斯-赛德尔迭代法8.1 算法说明步1 取初始点，精度要求，最大迭代次数N,置步2计算步3 假设那么停算，输出作为方程组的近似解；否那么转步 4.步4假设那么停算，输出迭代失败信息；否那么置，转步2.8.2 高斯-赛德尔迭代法流程图开

18、始 请输入精度要求和最大迭代次数 结束图8-1 高斯赛德尔迭代流程图8.3 高斯-赛德尔迭代法流运行结果图8-2 高斯赛德尔迭代运行结果9 复化辛普森公式求积分9.1 复化辛普森公式求积分算法说明将积分区间a,b分成n等分，分点为其中记区间的中点为，在每个小区间上用辛普森公式，那么可得到复化辛普森公式9.2 复化辛普森公式求积分流程图开始定义输入n,a,b输出结束图9-1 复化辛普森流程图9.3 复化辛普森公式运行结果图9-2 复化辛普森运行结果10 最小二乘线性拟合10.1 最小二乘线性拟合算法说明10-1设采集到一组N个样本数据，对这组数据进行线性拟合，求出拟合值。由每个样本数据做误差求和

19、计算：，err值由值确定，因此，误差和分别做偏微分计算，当每个偏微分值为0时，误差和最小。故有：10-210-3 根据每个偏微分等于零，可以得到两个关于的一次方程，求解该方程即可以分别得到的值。10.2 最小二乘法进行线性拟合流程图Start 输入自变量x与因变量y，数据点的坐标用x、y的对应次方构造方程解线性方程组得到a，bEnd图10-1 最小二乘进行线性拟合流程图10.3 最小二乘法进行线性拟合运行结果图10-2 最小二乘进行线性拟合运行结果11 设计体会及改良意见这次课程设计是我受益匪浅，掌握了很多的东西，不仅仅包括对于算法的进一步理解，还有对于c+程序编程能力的提高，以及对于工作软件

20、word的操作能力。在上课期间，由于理论知识比拟空泛、抽象，对于常用的算法没有理解透彻。所以，这次课程设计对我的帮助很大，通过编程更深一步加深对算法的理解与掌握，提高了我的上机实践能力。在实验的编程过程中遇到了很大的困惑，由于对算法理论的不熟悉以及C+学习的知识有限，编写的程序少且仅仅是一些小的算法。所以当数值算法与C+语言结合在一起的时候显得有些吃力。在最后的排版中也是非常的吃力，由于平时对于word操作不熟练导致这次课设进度缓慢。在这次课程设计编写的10个程序中存在一个问题：有的程序在解决数学问题时，是在程序中录入方程或者数据，而不是在运行界面输入。这样就会导致在更改问题之后，就必须在程序

21、代码中修改，以至于程序存在一定的局限性。这点是需要改良的。通过这次课程设计我深深地体会到数学方法和数学逻辑思维在计算机编程上的重要性，学好编程不仅仅是如何巧妙地实现问题，数学逻辑培养也是不容无视的。参 考 文 献1 谭浩强.C+程序设计.M.北京：清华大学出版社,2005.2 谭浩强.C+程序设计与上机指导(第三版).M.北京：清华大学出版社,2005.3数值方法MATLAB版马昌凤，林伟川编著.北京：科学出版社，2021.附 录附录一 高斯列主元法程序代码#include#include #includeusing namespace std;#define m 10/double Amm,

22、bm,Augmm+1; /全局数组void equ(double Amm+1,double xm,int n) int i,i1,j,k; double Augmm+1,maxele,Temp,l,s;for(i1=0;i1n-1;i1+)maxele=fabs(Augi1i1);k=i1;for(i=i1;in;i+) /找主元if(maxelefabs(Augii1)maxele=fabs(Augii1);k=i; cout检验找主元的正确性：主元位于第 k+1 行，第 i1+1 列 endlendl;for(j=i1;jn+1;j+) /实施逐个元素换行Temp=Augi1j;Augi1

23、j=Augkj;Augkj=Temp; /向用户输出提示信息 cout交换第 i1+1 行与第 k+1 行,请检验换行的正确性 Aug= endlendl; for(i=0;in;i+) /输出Aug，以供检验 for(j=0;jn+1;j+)coutsetw(10)Augij ; coutendl; coutendl;for(k=i1+1;kn;k+) /消元，以第i1行为工具行处理以下各行，下面代码描述了处理第k行的过程l=-Augki1/Augi1i1;for(j=i1;jn+1;j+) Augkj=Augkj+l*Augi1j;xn-1=Augn-1n/Augn-1n-1; /回带求解

24、cout在被调函数 equ 中，回带计算解X，边计算，边输出endlendl;cout xn=xn-1=0;i-)s=0;for(j=i+1;jn;j+)s=s+Augij*xj;xi=(Augin-s)/Augii; cout xi+1=xi; ; coutendl; /编写主函数int main() void equ(double Amm+1,double xm,int n);int i,j,n;static double Amm+1,bm,xm; cout请输入未知量的个数nn;if(nm)cout问题规模太大，需要修改源程序中等符常量m定义endlendl;return 0;for(i

25、=0;in;i+)cout请输入增广矩阵A的第i+1行:;for(j=0;jAij; /调用equ 函数解 n 元一次线性方程组 equ(A,x,n);/输出解Xcout在主函数中输出解x= ;coutsetiosflags(ios:fixed)setiosflags(ios:right)setprecision(4);for(i=0;in;i+)coutsetw(10)xi , ;coutendlendl;return 0;附录二 牛顿法解非线性方程组的程序代码#include#include#define N 2 / 非线性方程组中方程个数、未知量个数 #define Epsilon 0.

26、0001 / 差向量1范数的上限#define Max 100 /最大迭代次数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 x1

27、N); /由近似解向量 x0 计算近似解向量 x1float x0N,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,j,iter=0;cout初始近似解向量：endl;for (i=0;ix0i; coutendl;coutendl;doiter=iter+1;cout第 iter 次迭代开始endl; /计算向量函数的因变量向量 y0ff(x0,y0); /计算雅克比矩阵 jacobianffjacobian(x0,jacobian); /计算雅克比矩阵的逆矩阵 invjacobianinv_jacobian(jacobian,invjac

28、obian); /由近似解向量 x0 计算近似解向量 x1newdundiedai(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; int i; x=xx0; y=xx1; yy0=x*x-2*

29、x-y+0.5; yy1=x*x+4*y*y-4; /计算初值位置的值 cout向量函数的因变量向量是： endl; for( i=0;iN;i+) coutyyi ; coutendl; coutendl;void ffjacobian(float xxN,float yyNN) float x,y; int i,j; x=xx0;y=xx1;/jacobian have n*n element /计算函数雅克比的值yy00=2*x-2;yy01=-1;yy10=2*x;yy11=8*y;cout雅克比矩阵是： endl;for( i=0;iN;i+) for(j=0;jN;j+) cout

30、yyij ; coutendl; coutendl;void inv_jacobian(float yyNN,float invNN)float augNN2,L; int i,j,k; 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 (k=i+1;kN;k+) L=-augki/augii;for(j=i;j0;i-) for (k=i-1;k=0;k-) L=-augki/augii;for(j=N2-1;j=0;j-)

31、augkj=augkj+L*augij; for (i=N-1;i=0;i-) for(j=N2-1;j=0;j-)augij=augij/augii;for (i=0;iN;i+) for(j=N;jN2;j+) invij-N=augij;cout雅克比矩阵的逆矩阵： endl;for (i=0;iN;i+) for(j=0;jN;j+) coutinvij ; 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

32、+) sum=sum+invij*y0j;x1i=x0i-sum; coutendl; cout近似解向量：endl;for (i=0;iN;i+) coutx1i ; coutendl;coutendl;附录三 经典四阶龙格库塔法#include#includeusing namespace std;double f(double x,double y,double z); /函数声明double g(double x,double y,double z); /函数声明int main() double z0,x0,y0; /定义变量 couty0z0; int n; coutn; doub

33、le a,b,h; /求解的区间为0,5 a=0;b=5; h=(b-a)/n; cout输出步长:h=hendl; double f1,f2,f3,f4,g1,g2,g3,g4; double x1,z1,y1; x0=0; for(int k=0;kn;k+) /龙格库塔解一阶线性方程组f1=f(x0,y0,z0);g1=g(x0,y0,z0);f2=f(x0+h/2,y0+h/2*f1,z0+h/2*g1);g2=g(x0+h/2,y0+h/2*f1,z0+h/2*g1);f3=f(x0+h/2,y0+h/2*f2,z0+h/2*g2);g3=g(x0+h/2,y0+h/2*f2,z0+

34、h/2*g2);f4=f(x0+h,y0+h*f3,z0+h*g3);g4=g(x0+h,y0+h*f3,z0+h*g3);x1=x0+h;x0=x1; y1=y0+h/6*(f1+2*f2+2*f3+f4);y0=y1;z1=z0+h/6*(g1+2*g2+2*g3+g4);z0=z1; coutx0=x0t; couty0=y0t; coutz0=z0endl;return 0; double f(double x,double y,double z) /定义方程组第一个函数 double r; r=-0.01*y-99.99*z; return r; double g(double x,

35、double y,double z) /定义方程组第二个函数 double p; p=-100*z; return z; 附录四 三次样条插值算法压紧样条#include#includeusing namespace std;const int max = 50;float xmax, ymax, hmax;float cmax, amax, fxymmax;float f(int x1, int x2, int x3)float a=(yx3-yx2) / (xx3 - xx2);float b = (yx2 - yx1) / (xx2 - xx1);return (a-b)/(xx3-xx

36、1); /求差分void cal_m(int n) float Bmax;B0=c0 / 2;for(int i=1; in; i+)Bi=ci/(2-ai*Bi-1);fxym0=fxym0/2;for(i=1;i=0;i-)fxymi=fxymi-Bi*fxymi+1;void printout(int n);int main()int n,i; char ch;do coutn;for(i =0;i=n;i+)cout请输入Xixi; cout请输入Yiyi; for(i =0; i n; i+) /求步长hi=xi+1-xi;coutt;switch(t) case 1:cout输入 Y0 Ynf0f1;c0 = 1; an = 1;fxym0 = 6*(y1