数值计算课程设计

1、数值计算课程设计说明书题目： 典型数值算法的C+语言程序设计 院 （系）： 理学院 专业班级： 数学1xx 学 号： xxxx12010134 学生姓名： x 指导教师： xxxx 2014年7月 6日陕 西 科 技 大 学数值计算课程设计任务书理学院信息与计算科学/应用数学专业信息12*/数学12* 班级 学生： xxx 题目：典型数值算法的C+语言程序设计 课程设计从 2014 年 5 月 26 日起到 2014 年 7 月 6 日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，下2cm，左2.5cm、右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个图是关于高斯列主元法解方程

6、组算法图，图的下方的图号图名应为：图1-5 高斯列主元法解方程组算法图。（b）表格：表的名称及表内文字采用中文，中文字体为五号宋体，表名在表格上面。表号以章为单位顺序编号，表内必须按规定的符号标注单位。格式为：表1-1，空一字符后接表格名称。（c）公式：公式书写应在文中另起一行，居中排列。公式序号按章顺序编号。字体为五号宋体，序号靠页面右侧。格式为：（1-1）。设计体会及今后的改进意见：设计总结要写出算法理解，编程经验等技术性、学术性总结；体会要简洁、真实、深刻，切忌空话、大话，客套话和矫揉造作之词。改进意见要合理、中肯。参考文献的要求：另起一页，居中打印参考文献四字（四号黑体，段前段后1行）

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

8、组11.1算法说明11.2算例11.3程序代码12牛顿法解非线性方程组22.1算法说明22.1算例：32.2程序代码53经典四阶龙格库塔法解一阶微分方程组93.1算法说明93.2算例103.3程序代码104三次样条插值算法（压紧样条）124.1算法说明124.2算例124.3 C+程序代码134.4 matlab程序代码165龙贝格求积分算法175.1算法说明175.2算例175.3 程序代码176 M次多项式曲线拟合186.1算法说明186.2算例196.3程序代码197拉格朗日插值解多项式227.1算法说明227.2算例227.3程序代码238二分法求解非线性方程248.1算法说明248.

9、1算例248.2程序代码249 不动点迭代269.1算法说明269.2算例269.3程序代码2610复化梯形求积分公式2710.1算法说明2710.2算例2810.3程序代码2811设计体会29参 考 文 献30数值计算课程设计1高斯列主元法解线性方程组1.1算法说明将线性方程组做成增广矩阵，对增广矩阵进行行变换。对第1列元素，在第i行及以下的元素中选取绝对值最大的元素，将该元素最大的行与第i行交换，然后采用高斯消元法将新得到的消去第i行以下的元素。一次进行直到。从而得到上三角矩阵。再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。1.2算例课本99页例3.16，求解方程组的解运行结果如

10、图1-1所示。图1-11.3程序代码#include#includeusing namespace std;const N=20;float aNN;int m;int main()int i,j;int c,k,n,p,r;float xN,lNN,s,d;coutm;coutendl;cout请按顺序输入增广矩阵a:endl;for(i=0;im;i+)for(j=0;jaij;for(i=0;im;i+) for(j=i;jfabs(aii)?j:i; /*找列最大元素*/for(n=0;nm+1;n+) s=ain; ain=acn; acn=s; /*将列最大数防在对角线上*/for

11、(p=0;pm+1;p+)coutaipt;coutendl;for(k=i+1;km;k+) lki=aki/aii; for(r=i;r=0;i-) d=0;for(j=i+1;jm;j+)d=d+aij*xj;xi=(aim-d)/aii; /*求解*/cout该方程组的解为:endl;for(i=0;im;i+)coutxi=xit; /system(pause); return 0; 2牛顿法解非线性方程组2.1算法说明设已知，第一步计算函数第二步计算雅可比矩阵第三步求线性方程组的解。第四步计算下一点重复以上过程。2.1算例：设有非线性方程组（课本137页例3.32）设初始值，用牛顿

12、法计算。解：迭代3次后的近似解向量为1.9068,0.311219;具体运行结果如图2-1，2-2，2-3所示。 图2-1图2-2图2-12.2程序代码#include#include#define N 2 /非线性方程组中方程的个数、未知量的个数#define Epsilon 0.0001 /差向量的上限#define Max 100 /最大迭代次数using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN); /计算向量函数的因变量向量yyNvoid ffyakebi(float xxN,float

13、yyNN);/计算雅可比矩阵yyNNvoid inv_yakebi(float yyNN,float invNN);/计算雅可比矩阵的逆矩阵invvoid newton(float x0N,float invNN,float y0N,float x1N);/由近似解向量x0求近似解向量x1float x0N=2.0,0.25,y0N,yakebiNN,invyakebiNN,x1N,errornorm;int i,iter=0;cout初始近似解向量:endl;for(i=0;ix0i;for(i=0;iN;i+)coutx0i ;coutendl;coutendl;doiter=iter+1

14、;cout第iter次迭代开始endl;ff(x0,y0);ffyakebi(x0,yakebi);inv_yakebi(yakebi,invyakebi);newton(x0,invyakebi,y0,x1);errornorm=0;for(i=0;iN;i+)errornorm=errornorm+fabs(x1i);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=yy1;yy0=x*x-2

15、*x-y+0.5;yy1=x*x+4*y*y-4;cout向量函数的因变量向量是:endl;for(i=0;iN;i+)coutyyi ;coutendl;coutendl;void ffyakebi(float xxN,float yyNN)float x,y;int i,j;x=xx0;y=xx1;yy00=2*x-2;yy01=-1;yy10=2*x;yy11=8*y;cout雅可比矩阵是:endl;for(i=0;iN;i+)for(j=0;jN;j+)coutyyij ;coutendl;coutendl;void inv_yakebi(float yyNN,float invNN)

16、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;elseaugij=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=1;jN2;j+)augkj=augkj+L*augij;for(i=0;iN;i+)for(j=0;jN2;j+)co

17、utaugij ;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;f

18、or(i=0;iN;i+)for(j=0;jN;j+)coutinvij ;coutendl;coutendl;void newton(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.1算法说明4阶龙格-库塔方法(RK4)可模拟N=4的泰勒方法

19、的精度。这种算法可以描述为，自初始点开始，利用下面的计算方法生成近似序列3.2算例：用4阶龙格-库塔方法计算区间0.0,0.2上式（3）的数值解，采用10个子区间，步长h=0.02。（课本400页例9.15）求解运行结果如图3-1所示。初值图3-13.3程序代码#include#includeusing namespace std;void LG(double (*f)(double t,double x,double y),double (*g)(double t,double x,double y),double chuzhi3,double resu3,double h)double f

20、1,f2,f3,f4,g1,g2,g3,g4,t0,x0,y0,x1,y1;t0=chuzhi0;x0=chuzhi1;y0=chuzhi2;f1=f(t0,x0,y0);g1=g(t0,x0,y0);f2=f(t0+h/2,x0+h*f1/2,y0+h*g1/2);g2=g(t0+h/2,x0+h*f1/2,y0+h*g1/2);f3=f(t0+h/2,x0+h*f2/2,y0+h*g2/2);g3=g(t0+h/2,x0+h*f2/2,y0+h*g2/2);f4=f(t0+h,x0+h*f3,y0+h*g3);g4=f(t0+h,x0+h*f3,y0+h*g3);x1=x0+h*(f1+2

21、*f2+2*f3+f4)/6;y1=y0+h*(g1+2*g2+2*g3+g4)/6;resu0=t0+h;resu1=x1;resu2=y1;int main()double f(double t,double x,double y);double g(double t,double x,double y);double chuzhi3,resu3;double a,b,S;double t,step;int i;coutchuzhi0chuzhi1chuzhi2;coutab;coutstep;S=(b-a)/step;coutchuzhi0setw(10)chuzhi1setw(10)c

22、huzhi2endl;for(i=0;istep;i+)LG(f,g,chuzhi,resu,S);coutresu0setw(10)resu1setw(10)resu2endl;chuzhi0=resu0;chuzhi1=resu1;chuzhi2=resu2;return 0;double f(double t,double x,double y)double dx;dx=x+2*y;return(dx);double g(double t,double x,double y)double dy;dy=3*x+2*y;return(dy);4三次样条插值算法（压紧样条）4.1算法说明表5-

23、1 三次样条插值算法说明表策略描述包含和的方程(i)三次紧压样条，确定，（如果导数已知，这是“最佳选择”）(ii)natural三次样条（一条“松弛曲线”），(iii)外挂到端点(iv) 是靠近端点的常量，(v)在每个端点处指定，4.2算例求三次紧压样条曲线，以经过点（0,0），（1,0.5）,(2,2.0)和（3,1.5），且一阶导数的边界条件为S(0)=0.2和S(3)=-1。（课本222页例5.7）求解运行结果界面如下：图4-1 三次样条插值算法（压紧样条）Matlab作图如图4-2所示。 图4-24.3 C+程序代码#include#includeusing namespace std

24、;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-xx1); /求差分void cal_m(int n)/用追赶法求解出弯矩向量M float Bmax; B0 = c0 / 2; for(int i = 1; i n; i+) Bi = ci / (2 - ai*Bi-1); fxym0 =

25、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+) coutPlease put in Xixi; coutPlease put in 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*(

26、y1 - y0) / (x1 - x0) - f0) / h0; fxymn = 6*(f1 - (yn - yn-1) / (xn - xn-1) / hn-1; break; case 2:cout输入 Y0 Ynf0f1; c0 = an = 0; fxym0 = 2*f0; fxymn = 2*f1; break; default:cout不可用n;/待定 ; for(i = 1; i n; i+) fxymi = 6 * f(i-1, i, i+1); for(i = 1; i n; i+) ai = hi-1 / (hi + hi-1); ci = 1 - ai; an = hn-

27、1 / (hn-1 + hn); cal_m(n); coutn输出三次样条插值函数：n; printout(n); coutch; while(ch = y | ch = Y); return 0;void printout(int n)coutsetprecision(6); for(int i = 0; i n; i+)couti+1: xi , xi+1nt;coutSi+1 0) cout-t*(x - xi+1)3; else cout-t*(x - xi+1 0) cout + t*(x - xi)3; else cout - t*(x - xi)3; cout 0) cout-

28、 t*(x - xi+1); else cout- -t*(x - xi+1 0) cout + t*(x - xi); else cout - -t*(x - xi); coutendl; coutendl;4.4 matlab程序代码x1=0:0.01:1;y1=0.06*(x1 - 1).3 + 0.42*(x1 - 0).3 - 0.06*(x1 - 1) + 0.08*(x1 - 0);x2=1:0.01:2;y2=-0.42*(x2 - 2).3 - 0.62*(x2 - 1).3 - 0.08*(x2 - 2) + 2.62*(x2 - 1);x3=2:0.01:3;y3=0.6

29、2*(x3 - 3).3 + 0.06*(x3 - 2).3 - 2.62*(x3 - 3) + 1.44*(x3 - 2);X=0 1 2 3;Y=0 0.5 2 1.5;plot(x1,y1,x2,y2,x3,y3,X,Y,*)gtext(S1)gtext(S2)gtext(S3)5龙贝格求积分算法5.1算法说明生成的逼近表，并以为最终解来逼近积分逼近存在于一个特别的下三角矩阵中，第0列元素用基于个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为当时，程序在第行结束。5.2算例用龙贝格求积公式计算的值，运行结果如图5-1所示。图5-15.3 程序代码#includ

30、e#include#includeusing namespace std;#define f(x) (sin(x) /列举函数#define N 20 /区间等分数#define MAX 10 /最大迭代次数#define a 0 /所求积分的上下限#define b 1#define epsilon 0.0001 /精度double Romberg(double aa,double bb,long int n)int i;double sum,h=(bb-aa)/n;sum=0;for(i=1;in;i+)sum=sum+f(aa+i*h);sum=sum+(f(aa)+f(bb)/2;re

31、turn (h*sum);void main()int i;long int n=N,m=0;double T2MAX+1;T10=Romberg(a,b,n);n=n*2;for(m=1;mMAX;m+)for(i=0;i=m;i+)T0i=T1i;T10=Romberg(a,b,n);n=n*2;for(i=1;i=m;i+)T1i=T1i-1+(T1i-1-T0i-1)/(pow(2,2*m)-1);if(T0m-1-T1m)epsilon)coutT=T1mendl;return;6 M次多项式曲线拟合6.1算法说明设有N个点，横坐标是确定的。最小二乘抛物线的系数表示为（6-1）求解A

32、，B和C的线性方程组为6.2算例 根据4个点（-3，3），（0，1），（2，1）和（4，3），求解最小二乘抛物线。（课本211页例5.6）求解结果如图6-1所示。图6-16.3程序代码#include#include#define MAX 20using namespace std;int main()int n,m,i,j,k;void inv(double XMAXMAX,int n,double EMAXMAX);double XMAX=0,YMAX=0,ZMAXMAX=0,BMAX=0,CMAX=0;double AMAXMAX=0,BZMAXMAX=0,EMAXMAX=0;coutn

33、;coutendl;cout请输入n个点的X坐标序列:;for(i=0;iXi;coutendl;cout请输入n个点的Y坐标序列:;for(i=0;iYi;coutendl;coutm;for(i=0;in;i+)for(k=1;k=m+1;k+)Zik-1=pow(Xi,k-1);/求Z的转置for(i=0;in;i+)for(j=0;jm+1;j+)BZji=Zij; /计算其转置的BF与Z的乘 for(i=0;im+1;i+)for(j=0;jm+1;j+)for(k=0;kn;k+)Aij+=BZik*Zkj;/计算Z的转置BZ与Y的乘for(i=0;im+1;i+)for(j=0;

34、jn;j+)Bi+=BZij*Yj;/调用inv函数求解矩阵A的逆矩阵Einv(A,n,E);/计算A的逆BZ与B的乘for(i=0;im+1;i+)for(j=0;jn;j+)Ci+=Eij*Bj;coutendl;cout拟合后的m次多项式系数为，幂次由高到低:=0;i-)coutCit;coutendl;cout拟合后的m次多项式为:endl;cout=0;i-)if(i=0)cout+Ci;elsecout+Ci*xi;coutendl;return 0;void inv(double XMAXMAX,int n,double EMAXMAX)/求解任意可逆矩阵的逆,X为待求解矩阵,E

35、为全零矩阵,非单位矩阵,也可以是单位矩阵int i,j,k;double temp=0;for(i=0;iMAX;i+)for(j=0;jMAX;j+)if(i=j)Eij=1;for(i=0;in-1;i+)temp=Xii;for(j=0;jn;j+)Xij=Xij/temp;Eij=Eij/temp;for(k=0;kn;k+)if(k=i)continue;temp=-Xii*Xki;for(j=0;jn;j+)Xkj=Xkj+temp*Xij;Ekj=Ekj+temp*Eij;7拉格朗日插值解多项式7.1算法说明7.2算例给定函数四个点坐标(1.1,3.887),(2.3,4.276

36、),(3.9,4.651),(5.1,2.117),试用拉格朗日插值确定函数在x=2.101处的值。求解运行结果界面如图7-1所示，可知函数在点2.101处的值为4.14569。图7-17.3程序代码#include#define N 4/插值节点的个数using namespace std;void main()float xN,yN,a,f=0,temp;int i,j;cout输入插值节点的坐标：endl;for(i=0;ixi;cinyi;couta;for(i=0;iN;i+)temp=1;for(j=0;jN;j+)if(i!=j)temp=temp*(a-xj)/(xi-xj);

37、f=f+temp*yi;cout所求值为：fendl;8二分法求解非线性方程8.1算法说明若要求已知函数的根(x的解),则第一步先找出一个区间a,b,使得与异号,因此这个区间内一定包含着方程的根。第二步求该区间的中点 ，并找出的值。若与正负号相同则取m,b为新的区间, 否则取a,m,重复第二和第三步至理想精确度为止。8.1算例利用二分法寻找函数在区间0，2内的的零点。（课本44页例2.7）运行结果如图8-1所示。图8-18.2程序代码#include #include #define eps 0.001double fun(double x)return x*sin(x)-1;double d

38、ichotomy(double a,double b)double c=0.0;if(fun(a)0)while(true)c=(a+b)/2;if(fun(c)0)a=c;if(fabs(a-b)eps)return (a+b)/2;else if(fun(c)=0)return c;elseb=c;if(fabs(a-b)eps)return (a+b)/2;elsecout你输入的a和b不正确endl;return -1;int main()double a=0;double b=2;double result=dichotomy(a,b);cout求解的结果是:resultendl;return 0;9 不动点迭代9.1算法说明 先将改写成，然后对进行迭代，即 其中然后判断是否成立，成立则返回，不成立就重复以上步骤9.2算例当使用函数时，设迭代。通过求解方程可找到不动点。(课本36页例2.4)求解运行结果界面如图9-1所示，可知迭代在区间-2，2上最终收敛到2。图9-19.3程序代码#include#include#include#define MAX 20#define eps 0.00001using