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

数值计算课程设计说明书题目： 典型数值算法的C+语言程序设计 院 （系）： 理学院 专业班级： 数学101班 学 号： 201012010112 学生姓名： 仰春兰 指导教师： 刘海峰 2012 年 6 月 21 日陕 西 科 技 大 学理学院信息与计算科学/应用数学专业信息10/数学101 班级 学生： 仰春兰 题目：典型数值算法的C+语言程序设计 课程设计从 2012 年 5 月 21 日起到 2012 年 6 月 27 日1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）：每人需作10个算法的程序、必做6题、自选4题。对每个算法要求用C+语言进行编程。必选题：1、经典四阶龙格库塔法解一阶微分方程组2、高斯列主元法解线性方程组3、牛顿法解非线性方程组4、龙贝格求积分算法5、三次样条插值算法（压紧样条）用C+语言进行编程计算 依据计算结果，用Matlab画图并观察三次样条插值效果。6、M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果。自选题：自选4道其他数值算法题目.每道题目重选次数不得超过5次.2、对课程设计成果的要求包括图表、实物等硬件要求：1）提交课程设计报告1. 按照算法要求，用C+语言设计和开发应用程序，提交由算法说明；程序设计说明（程序中的主要变量和数据结构说明）；系统技术文档 (包括系统各模块主要流程图，软件测试方案与测试记录、软件调试和修改记录、测试结论、运行情况记录)，系统使用说明书，源程序代码为附录构成的课程设计报告。2）课程设计报告版式要求目录的要求：居中打印目录二字，（四号黑体，段后1行），字间空一字符；章、节、小节及其开始页码（字体均为小四号宋体）。节向右缩进两个字符，小节及以后标题均向右缩进四个字符。目录中应包含正文及其后面部分的条目。目录的最后一项是无序号的“参考文献资料”正文的要求：计算正确，论述清楚，文字简练通顺，插图简明，书写整洁。文中图、表按制图要求绘制。打印版面要求：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）顺序表示，如果编号及其后内容新起一个段落，则编号前空两个中文字符。曲线图表要求：所有曲线、图表、线路图、流程图、程序框图、示意图等不准徒手画，必须按国家规定标准或工程要求绘制（应尽可能采用计算机辅助绘图）。课程设计说明书（报告）中图表、公式要求如下：（a）图：图的名称采用中文，中文字体为五号宋体，图名在图片下面。引用图应在图题右上角标出文献来源。图号以章为单位顺序编号。格式为：图1-1，空一字符后接图名。（b）表格：表的名称及表内文字采用中文，中文字体为五号宋体，表名在表格上面。表号以章为单位顺序编号，表内必须按规定的符号标注单位。格式为：表1-1，空一字符后接表格名称。（c）公式：公式书写应在文中另起一行，居中排列。公式序号按章顺序编号。字体为五号宋体，序号靠页面右侧。格式为：（1-1）。设计体会及今后的改进意见：体会要简洁、真实、深刻；忌空话、套话等不实之词。改进意见要合理、中肯。参考文献的要求：另起一页，居中打印参考文献四字（四号黑体，段前段后1行），字间空一字符；另起一行，按论文中参考文献出现的先后顺序用阿拉伯数字连续编号（参考文献应在正文中注出）；参考文献中每条项目应齐全（字体均为小四号宋体）。（格式：编号作者.论文或著作名称.期刊名或出版社.出版时间）。（期刊应注明第几期、起止页数（包括论著）。3）设计报告装订顺序与规范封面数值计算课程设计任务书目录数值计算设计课程设计报告正文设计体会及今后的改进意见参考文献（资料）左边缘装订3、课程设计工作进度计划：时间设计任务及要求第14周编写和调试程序并按要求撰写设计报告 指导教师： 日期： 教研室主任： 日期： 目 录第1节 高斯列主元法解线性方程组1.1 高斯列主元法概要 首先，写出非线性方程组的增光矩阵； 其次，对增广矩阵进行变换；其中包括；一，找出每一列的最大元，并与对角线元素进行兑换；二，进行高斯列主元消元，用对角线元素将其所在列下其他元素变为0，直到得到一个上三角矩阵；最后，对上述所得矩阵进行回代解线性方程组，即可以得到方程组的解。1.2 算法流程图1.3 程序调试输入初值:4 0.01 0.01p=1k=2err=3relerr=2.9703p=1.75k=3err=0.75relerr=0.426136 p=1.98438k=4err=0.234375relerr=0.117518p=1.99994k=5err=0.015564relerr=0.0077435p=1.99994k=5err=0.015564Press any key to continue- 8 -1.4 程序代码#include#define n 4/代表有4个未知量n个方程using namespace std;int main() int i,j,t,p,g,r;double l,ann+1,max,Q,xn; cout输入增广矩阵：endl;/输入增广矩阵 for(i=0;in;i+) for(j=0;jaij; coutnow array a: endl; for(i=0;in;i+) for(j=0;jn+1;j+) coutaij ; coutendl; for(p=0;pn;p+) max=abs(app);r=p;/假设对角线元素为最大 for(i=p+1;imax)/将第i行p列的元素与主对角线的元素进行比较for(j=p;j5;j+)g=apj;apj=aij;aij=g; r=i; /将主对角线元素所在行与最大元素所在行进行对换 /高斯列主元消元l=aip/app; for(t=p;tn;t+) ait=ait-l*apt; coutnow array a: endl; for(i=0;i4;i+) for(j=0;j5;j+) coutaij ; coutendl; /输出现在的增广矩阵 cout回代:endl;xn-1=an-1n/an-1n-1;coutxn-1=xn-1=0;i-) Q=0;for(j=i+1;jn;j+) Q=Q+xj*aij;xi=(ain-Q)/aii;coutxi=xi=0) return z; else return -z;/求绝对值 第2节 经典四阶龙格库塔法阶一阶微分方程组2.1 算法概要可模拟N=4的泰勒方法的精度，这种算法可以描述为，自初始点 开始，利用生成的近似序列，其中2.2 算法流程开始输入x0,y0,a,b,nh=(b-a)/2k=0f1=f(t0,x0,y0)g1=g(t0,x0,y0)f2=f(t0+h/2,x0+h/2*f1,y0+h/2*g1)g2=g(t0+h/2,x0+h/2*f1,y0+h/2*g1)f3=f(t0+h/2,x0+h/2*f2,y0+h/2*g2)g3=g(t0+h/2,x0+h/2*f2,y0+h/2*g2)f4=f(t0+h,x0+h*f3,y0+h*g3)g4=g(t0+h,x0+h*f3,y0+h*g3)输出x0,y0,t0,iink=k+1结束2.3 程序调试输入初值:6 4输入迭代次数:100输入步长:h=0.05t0=-9.25596e+061x0=1.94022e+009y0=2.91033e+009Press any key to continue2.4程序代码#include#includeusing namespace std;float f(double t,double x,double y);float g(double t,double x,double y);int main() double t0,x0,y0; coutx0y0; int n; coutn; double a,b,h; a=0;b=5; h=(b-a)/n; cout输入步长:h=hendl; double f1,f2,f3,f4,g1,g2,g3,g4; double x1,t1,y1; for(int k=0;kn;k+)f1=f(t0,x0,y0);g1=g(t0,x0,y0);f2=f(t0+h/2,x0+h/2*f1,y0+h/2*g1);g2=g(t0+h/2,x0+h/2*f1,y0+h/2*g1);f3=f(t0+h/2,x0+h/2*f2,y0+h/2*g2);g3=g(t0+h/2,x0+h/2*f2,y0+h/2*g2);f4=f(t0+h,x0+h*f3,y0+h*g3);g4=g(t0+h,x0+h*f3,y0+h*g3);t1=t0+h;t0=t1; x1=x0+h/6*(f1+2*f2+2*f3+f4);x0=x1;y1=y0+h/6*(g1+2*g2+2*g3+g4);y0=y1; coutt0=t0 endl; coutx0=x0 endl; couty0=y0 K的逼近表，比并以R(J+1，J+1)为最终解来逼近积分,在区间0,1上的积分，逼近存在于一个特别的下三角矩阵中，第0列元素用基于个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为当时，程序在第行结束。3.2算法流程图3.3程序调试输入初值:1 2 0.1输出矩阵:2.5 0 0 02.375 2.33333 0 02.34375 2.33333 2.33333 02.33594 2.33333 2.33333 2.33333h=0.0625quad=2.33333err=0Press any key to continue3.4程序代码#include#includeusing namespace std;float f(double x);int main()double a,b,delta,h,quad,err,s,x;double R55=0;int n,i;coutabndelta;int m,k,j,p;m=1;h=b-a;err=1;j=0;R00=h*(f(a)+f(b)/2;while(errdelta)&(jn)|(j4) s=0;j=j+1; h=h/2;for(int p=1;p=m;p+)x=a+h*(2*p-1);s=s+f(x);Rj0=Rj-10/2+h*s;m=2*m; for(int k=0;kj;k+)Rjk+1=Rjk+(Rjk-Rj-1k)/(pow(4.0,k+1)-1); err=abs(Rj-1j-1-Rjk+1);quad=Rjj;cout输出矩阵:;for(i=0;i4;i+)for(j=0;j4;j+)coutRij ;coutendl;couth=hendl;coutquad=quadendl;couterr=errendl;return 0; float f(double x) return x*x; 第4 节 牛顿法解非线性方程组4.1 算法概要设已知。第1步：计算函数第2步：计算雅可比矩阵第3步：求线性方程组的解。第4步：计算下一点重复上述过程。42算法流程图4.3 程序调试输入增广矩阵：1 2 1 4 132 0 4 3 284 2 2 1 20-3 1 3 2 6输出矩阵:4 2 2 1 200 -1 3 2.5 180 1.5 0.5 3.75 80 2.5 4.5 2.75 21输出矩阵:4 2 2 1 200 2.5 4.5 2.75 210 0 -2.2 2.1 -4.60 0 4.8 3.6 26.4输出矩阵:4 2 2 1 200 2.5 4.5 2.75 210 0 4.8 3.6 26.40 0 0 3.75 7.5输出上三角矩阵:4 2 2 1 200 2.5 4.5 2.75 210 0 4.8 3.6 26.40 0 0 3.75 7.5回代:X3=2X2=4X1=-1X0=3Press any key to continue4.4程序代码#include#includeusing namespace std;#define n 4int main()int i,j,p,k ,t;double m,Augnn+1,Xn,temp,max,s;cout输入增广矩阵：endl;for(i=0;in;i+) for(j=0;jAugij;/*for(i=0;i4;i+) for(j=0;j5;j+) coutAugij;*/for(p=0;pn-1;p+) max=fabs(Augpp);t=p;for(i=p+1;imax) max=fabs(Augip); t=i; for(j=0;jn+1;j+) temp=Augpj; Augpj=Augtj; Augtj=temp; for(k=p+1;kn;k+) m=Augkp/Augpp; for(j=p;jn+1;j+)Augkj=Augkj-m*Augpj;cout输出矩阵:endl; for(i=0;in;i+) for(j=0;jn+1;j+) cout Augij ; coutendl; cout输出上三角矩阵:endl; for(i=0;in;i+)for(j=0;jn+1;j+) cout Augij ; coutendl;cout回代:endl;Xn-1=Augn-1n/Augn-1n-1;coutXn-1=Xn-1=0;i-) s=0;for(j=i+1;jn;j+) s=s+Xj*Augij;Xi=(Augin-s)/Augii;coutXi=Xiendl; return 0;第5节 三次样条插值算法5.1 算法概要策略描述包含和的方程(i)三次紧压样条，确定，（如果导数已知，这是“最佳选择”）(ii)natural三次样条（一条“松弛曲线”），(iii)外挂到端点(iv) 是靠近端点的常量，(v)在每个端点处指定，5.2 程序调试求解经过点（1，0.0），（2，0.5），（3，2.0）和（4，1.5）而且一阶导数边界条件S(0)=0.2和S(3)=-1的三次压紧样条曲线求得的三次压紧样条曲线的矩阵S为: 0.0928571 0.45 -0.0428571 0 -1.02857 0.728571 1.8 0.5 1.92857 -2.35714 -0.0714286 2在区间（0,1)上的样条为：y=+0.0928571x3+0.45x2+-0.0428571x+0在区间（0,1)上的样条为：y=+-1.02857x3+0.728571x2+1.8x+0.5在区间（0,1)上的样条为：y=+1.92857x3+-2.35714x2+-0.0714286x+25.3 程序代码#include#include#define MAX 4double *diff(double X,int n) int i=0;double *H=NULL;H=(double*)malloc(n-1)*sizeof(double);for(i=1;i=n-1;i+)Hi-1=Xi-Xi-1;return H;double *divide(double Y,int N,double H)int i=0;double *D=NULL;D=(double*)malloc(N*sizeof(double);for(i=0;iN;i+)Di=Yi/Hi;return D;main()double XMAX=0,1,2,3,YMAX=0,0.5,2.0,1.5,SMAXMAX=0,temp=0.0,MMAX=0;int N=MAX-1,i=0,k=0;double AMAX-1-2=0,BMAX-1-1=0,CMAX-1-1=0;double dx0=0.2,dxn=1.0;double *H=NULL,*D=NULL,*U=NULL;printf(求解经过点（0，0.0），（1，0.5），（2，2.0）和（3，1.5）n而且一阶导数边界条件S(0)=0.2和S(3)=-1的三次压紧样条曲线nn);H=diff(X,MAX);D=divide(diff(Y,MAX),N,H);for(i=1;iN-2;i+)Ai=Hi+1;for(i=0;iN-1;i+)Bi=2*(Hi+Hi+1);for(i=1;iN-1;i+)Ci=Hi+1;U=diff(D,N);for(i=0;iN;i+)Ui=Ui*6;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(k=2;k=1;k-) Mk=(Uk-1-Ck-1*Mk+1)/Bk-1M0=3*(D0-dx0)/H0-M0/2;MN=3*(dxn-DN-1)/HN-1-MN-1/2;for(k=0;k=N-1;k+)Sk0=(Mk+1-Mk)/(6*Hk);Sk1=Mk/2;Sk2=Dk-Hk*(2*Mk+Mk+1)/6;Sk3=Yk;printf(求得的三次压紧样条曲线的矩阵S为:n);for(i=0;iMAX-1;i+)for(k=0;kMAX;k+)printf(%lft,Sik);printf(n);for(i=0;iN;i+)printf(n在区间(0,1)上的样条为:y=%+lfx3%+lfx2%+lfx%+lf,Si0,Si1,Si2,Si3);getch();第6节 M次多项式曲线拟合6.1 算法概要设有N个点，横坐标是确定的。最小二乘抛物线的系数表示为求解A，B和C的线性方程组为6.2 程序调试输入X:0.25 1.0 1.5 2.0 2.4 5.0输入Y:23.1 1.68 1.0 0.84 0.826 1.2576输出矩阵:1 0.25 0.06251 1 11 1.5 2.251 2 41 2.4 5.761 5 25输出矩阵:1 1 1 1 1 10.25 1 1.5 2 2.4 50.0625 1 2.25 4 5.76 25输出矩阵:6 12.15 38.072512.15 38.0725 151.21538.0725 151.215 680.244shu chu ju zhen:6 12.15 38.0725 1 0 012.15 38.0725 151.215 0 1 038.0725 151.215 680.244 0 0 1 6 12.15 38.0725 1 0 04.76837e-007 13.4688 74.1178 -2.025 1 00 -2.87507e-006 30.7918 4.79805 -5.50294 16 3.58696e-007 -3.17046e-006 7.31257 -6.04697 0.9349324.76837e-007 13.4688 3.48805e-006 -13.5742 14.2459 -2.407060 -2.87507e-006 30.7918 4.79805 -5.50294 11 5.97827e-008 -5.2841e-007 1.21876 -1.00783 0.1558224.76837e-007 1 2.58973e-007 -1.00783 1.0577 -0.1787140 -2.87507e-006 1 0.155822 -0.178714 0.0324762A 的逆矩阵1.21876 -1.00783 0.155822-1.00783 1.0577 -0.1787140.155822 -0.178714 0.0324762C0=22.9308C1=-16.9619C2=2.55319Press any key to continue6.3程序代码#include#include#define M 2#define n 6using namespace std;const int N2=2*(M+1);int main()float Xn,FnM+1,AM+1M+1,Yn,BM+1,CM+1,F0M+1n,ansM+1M+1;int i,j,k,p;float ss,s,m;cout输入X:;for(i=0;iXi;/gei ju zhen F fu zhi cout输入Y:;for(i=0;iYi; for(k=0;kM+1;k+) for(i=0;in;i+)Fik=pow(Xi,k);cout输出矩阵:endl; for(i=0;in;i+) for(j=0;jM+1;j+) coutFij ; coutendl; coutendl; /求矩阵F的逆矩阵F0for(i=0;iM+1;i+)for(j=0;jn;j+) F0ij=Fji;/输出F的逆矩阵 cout输出矩阵:endl; for(i=0;iM+1;i+) for(j=0;jn;j+) coutF0ij ; coutendl; coutendl; /求矩阵F与其转置的乘积，并将其赋值给矩阵A for(i=0;iM+1;i+)for(p=0;pM+1;p+)ss=0; for(j=0;jn;j+) ss=ss+F0ij*Fjp; Aip=ss; cout输出矩阵:endl; for(i=0;iM+1;i+) for(j=0;jM+1;j+) coutAij ; coutendl; coutendl;/qiu ju zhen B for(i=0;iM+1;i+) s=0; for(j=0;jn;j+) s=s+(F0ij*Yj); Bi=s; /求矩阵A的逆矩阵 float augM+1N2;for(i=0;iM+1;i+)for(j=0;jM+1;j+)augij=Aij;for(i=0;iM+1;i+)for(j=M+1;jN2;j+) augij=0; augiM+1+i=1;coutshu chu ju zhen:endl; for(i=0;iM+1;i+)for(j=0;jN2;j+)coutaugij ; coutendl;coutendl;/hua xia san jiao wei 0,bu han dui jiao xianfor(p=0;pM+1;p+) for(i=p+1;iM+1;i+) m=-augip/augpp; for(j=p;jN2;j+) augij=augij+m*augpj; /shu chu augfor(i=0;iM+1;i+)for(j=0;jN2;j+)coutaugij ;coutendl;cout=0;i-) for(k=i-1;k=0;k-)m=-augki/augii;for(j=N2-1;j=0;j-)augkj=augkj+m*augij; /shu chu augfor(i=0;iM+1;i+)for(j=0;jN2;j+)coutaugij ;coutendl;cout=0;i-) for(j=N2-1;j=0;j-) augij=augij/augii; /shu chu augfor(i=0;iM+1;i+)for(j=0;jN2;j+)coutaugij ;coutendl; for(i=0;iM+1;i+) for(j=M+1;jN2;j+)ansij-M-1=augij; coutA 的逆矩阵 endl; for (i=0;iM+1;i+) for(j=0;jM+1;j+) coutansij ; coutendl; coutendl;coutendl;for(i=0;iM+1;i+) double ss=0; for(j=0;jM+1;j+) ss=ss+ ansij*Bj; Ci=ss; coutCi=Ci00输入a，b输入a，bC=(a+b)/2;yc-f(c);yc=0ya*yc0b=ca=cx=c输出x的值 结束7.3程序调试输入初值:-2 2 0.0131233123无解a=-2b=2c=-9.25596e+061err=-1.07374e+008yc=-1.07374e+008Press any key to continue7.4程序代码#include#includeusing namespace std;double f(double x);int main()float err; int i,max1; double a,b,c,delta; float ya,yb,yc;coutabdelta; ya=f(a);yb=f(b);if(f(a)*f(b)0) cout无解endl;elsemax1=1+(int)(log(b-a)-log(delta)/log(2);/最大迭代次数for(i=1;i=max1;i+)c=(a+b)/2;yc=f(c);if(yc=0) a=c;b=c;coutc=c endl; coutyc=yc 0) b=c;yb=yc; else a=c;ya=yc;err=fabs(b-a); coutb=b endl; coutc=c delta) break;couta=a endl;coutb=b endl;coutc=c endl;couterr=err endl;coutyc=yc endl;return 0;double f(double x) double yx; yx=x*x-2*x+23;coutyx endl;return yx;第8节 牛顿-拉夫森解非线性方程8.1 算法概要利用求解非线性方程。8.2 程序调试输入初值:-1.5 0.1 0.01 50err=0.65relerr=13fabs(y)=0.4225err=0.114189relerr=2.28378fabs(y)=0.0130392err=0.00375594relerr=0.0751187-fabs(y)=1.41392e-005err=0.00375594k=2y=1.41392e-005p0=-0.732055Press any key to continue8.3程序代码#include#includeusing namespace std;int main() float ff(float x); float dff(float x); double err,relerr,delta,epsilon; int max1,k; float p0,p1,y; coutp0deltaepsilonmax1; for(k=0;k=max1;k+) p1=p0-ff(p0)/dff(p0); err=fabs(p1-p0); relerr=2*err/(fabs(p1)+delta); p0=p1; y=ff(p0); couterr=err endl; coutrelerr=relerr endl; coutfabs(y)=fabs(y) endl; if(errdelta)|(relerrdelta)|(fabs(y)epsilon) break; couterr=err endl; coutk=k endl; couty=y endl; coutp0=p0 endl; return 0;float ff(float x) float f(x); return x*x-2*x-2;float dff(float x) float df(x); return 2*x-2;第9