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

1、 目 录1用连续梯形公式求解积分11.1算法理解11.2具体的算例12高斯选主元解线性方程组22.1算法理解22.2算例及运行界面23经典四阶龙格库塔法解一阶微分方程组33.1算法理解33.2算例及运行界面34三次样条插值算法（紧压样条）54.1算法理解54.2具体算例及运行结果55龙贝格求积分算法75.1算法理解75.2算例及运行结果76 M次多项式曲线拟合86.1算法分析86.2算例及运行结果87牛顿 拉夫森迭代法97.1算法分析97.2算例及运行结果98牛顿法解非线性方程组108.1算法分析108.2算例及结果分析109 雅可比迭代解方程129.1算法分析129.2算例及结果分析1210

2、 指数拟合1410.1算法分析1410.3算例及结果显示1411 体会及总结15参考文献16附 录17数值计算课程设计1用连续梯形公式求解积分1.1 算法理解连续梯形公式可以求解定积分，其中a,b是积分函数的上下限，对于不同的积分，只需要变动被积函数和积分的上下限即可。数据类型使用双精度类型得到的结果更加准确。此程序通过利用指向函数的指针，可以计算各种函数的定积分，其中也包括不能求得原积分函数的积分，具有普遍的积分意义，若要对不同的积分函数进行积分的话，只需变动被积函数即可。1.2具体的算例 被积函数=x3;上下限分别为0,1. 输入参数，得到运行结果如下：图1-1 运行界面被积函数在0,1上

3、的利用梯形公式求解积分的运行界面：图1-2 运行界面37数值计算课程设计2高斯选主元解线性方程组2.1算法理解高斯选主元法可以很好的避免计算中分母为零的情况，虽然算法较一般的高斯消元法步骤更多，但程序相对而言更加完善，程序运行时可以更大程度的避免出错，若线性方程组出现无限多解的情况时，可以通过选主元的方法极早发现并给出错误提示。其中最主要的思想就是高等代数里面的矩阵的等价的初等变换。逐步的将矩阵化简成上三角型阵，最后从下至上，从后到前，进行一步步的回代求解。直到完全求解出方程组的解为止。2.2算例及运行界面例如可以利用编写好的程序解线性方程组 (2-1)运行结果如下所示：图2-1 高斯选主元运

4、行界面通过函数调用求解得到线性方程组的解x=3 -1 4 2。数值计算课程设计3经典四阶龙格库塔法解一阶微分方程组3.1 算法理解通过，可以得： （3-1） （3-2）迭代的思想与高斯赛德尔迭代类似，都充分的利用到了刚计算出的数据，数据的利用率高。，3.2算例及运行界面定义函数，。二元微分方程组的一般形式为：，将初值带入后。可以得到x,y关于参数t的等式，例如本例的迭代式为，通过输入初值并进行迭代，可以得到运行结果如下：图3-1 运行界面数值计算课程设计4三次样条插值算法（紧压样条）4.1算法理解三次样条插值通过插值点的连续性、一阶可导性以及二阶可导性和端点等约束来构造等式条件，通过求解三角线

5、性方程组H=MV，解得,最后解出每两个插值点之间的插值函数的三次多项式的系数，还可以用插值函数去预测未知点处的函数值。 4.2具体算例及运行结果得到运行结果如下：图4-1运行界面依据计算结果，用Matlab画图并观察三次样条插值效果，并且可以利用求解出的插值函数去预测某一点处的值。通过编写matlab程序并运行可以得： X=0 1 2 3; Y=0 0.5 2 1.5; dx0=0.2; dxn=-1; S=csfit(X,Y,dx0,dxn)S = 0.4800 -0.1800 0.2000 0 -1.0400 1.2600 1.2800 0.5000 0.6800 -1.8600 0.68

6、00 2.0000 x1=0:0.01:1; y1=polyval(S(1,:),x1-X(1); x2=1:0.01:2; y2=polyval(S(2,:),x2-X(2); x3=2:0.01:3; y3=polyval(S(3,:),x3-X(3); plot(x1,y1,x2,y2,x3,y3,X,Y,.) plot(x1,y1,x2,y2,x3,y3,X,Y,.) y=polyval(S(2,:),1.5-X(2)y = 1.3250图4-2 三次样条拟合曲线 数值计算课程设计1. 5龙贝格求积分算法5.1算法理解龙贝格积分通过不断缩短步长，是梯形积分，辛普森积分以及布尔积分的改进

7、。数据类型采用双精度度型，精度更高，也可以直接定义被积函数，然后利用函数指针实现调用。5.2算例及运行结果 可以求得积分的近似值。以下是函数在1,2的龙贝格积分所计算出的值。其中运行的最终结果如下：图5-1 在1,2上的积分值以上求解的界面是cos(x)的积分,当然也可以在一定范围内通过变换参数及函数求解其他函数的定积分。 数值计算课程设计6 M次多项式曲线拟合6.1算法分析 多项式拟合是通过最小二乘来确定多项式的系数，根据已知点的个数以及坐标，可得，通过系数求偏导的值为零，可以得到一个关于系数a,b,c的线性方程组，再利用高斯消元法或者高斯选主元解出多项式的系数。6.2算例及运行结果输入的四

8、个点分别是-3,3,0,1,2,1,4,3,其中拟合多项式的次数为2次，运行结果如下：图6-1 拟合界面通过matlab画图可以得到拟合图像与点之间的关系： X=-3 0 2 4; Y=3 1 1 3; x=-5:0.01:5; y=0.850519-0.192495*x+0.178462*x.2; plot(X,Y,+,x,y)图6-2 拟合曲线 数值计算课程设计7牛顿 拉夫森迭代法7.1算法分析牛顿-拉夫森法通过初始点处的切线与x轴的交点，不断迭代使得点趋近方程的解，其中迭代公式为，也即，迭代的要求是f(x)的一阶及其二阶函数都必需连续。通过这种方法可以求得部分非线性方程的解，比如说可以求

9、得一些无理数的近似值等。例如。7.2算例及运行结果其中求解的函数为，初始迭代点为x0=1.5,得到程序的运行结果为：图7-1 牛顿-拉夫森法求解方程根的迭代界面 牛顿法解非线性方程组8牛顿法解非线性方程组8.1算法分析 ,设，通过分别对u,v求偏导可以得到雅克比矩阵的具体形式，然后利用，解得，再利用进行迭代。对于多个未知数的方程组，则程序中需要变动函数的形式以及雅克比矩阵的形式。同时变换n的值。8.2算例及结果分析非线性方程组为： （8-1）运行情况如下：初始点为（1.2,1.2）图8-1 运行结果界面换初始迭代点后的运行情况：（0.2,0.2）图8-2 运行结果界面因此非线性方程组的迭代结果

10、是与初始点密切相关的，因为解的个数不唯一，不确定到底会收敛到哪一个具体的点。 雅克比迭代解方程9 雅可比迭代解方程9.1算法分析 线性方程组除了可以用高斯迭代法求解之外，还可以用其他的方法求解，比如说雅克比迭代，其中每一个方程都可以通过变化得到一个迭代式，通过给定一个初始点，可以进行迭代求得方程组的解，但是，雅克比迭代不一定收敛，当方程组的系数矩阵不严格对角占优时，利用雅克比迭代不一定能够得到方程组的解，交换方程组中方程的顺序都会影响迭代的结果。9.2算例及结果分析线性方程组为： （9-1）当方程组的系数矩阵不严格对角占优时，利用雅克比迭代不一定能够得到方程组的解，交换方程组中方程的顺序都会影

11、响迭代的结果。图9-1 运行结果界面图9-2 运行结果界面 迭代次数超过上限，因此雅克比迭代不一定都收敛。 数值计算课程设计10 指数拟合10.1算法分析指数拟合，即对于给定的一系列已知点，用一个指数函数去拟合，最后求解出具体的指数函数，与一般的拟合有相似之处，以前是用M次的多项式去拟合已知的点，而现在只是用指数类型的函数去拟合已知点而已。10.3算例及结果显示已知的需要拟合的点分别是1,2，2,6，3,9，通过指数拟合，可以得到一个确切的指数函数。程序的运行运行界面及结果如下：图10-1 运行结果界面指数拟合，即用一个指数函数y=a*exp(b*x)去拟合已知的点11 体会及总结 这次课程设

12、计主要的内容是将这学期所学的基本算法中的十个matlab程序再次编写成C+程序，并在Visual C+ 6.0上进行调试，修改和运行。通过老师的指导以及一周的练习，我深刻地意识到不同学科之间的联系原来是如此的紧密，最主要的是要学会如何思考和解决问题，学会怎么去分析和解决实际问题，比如说怎样实现程序在不同语言间的转换，这不仅要考虑到不同语言的特点和不同，更需要通过掌握不同语言的语法特性，进而实现转换。例如C+中的数组的下标是从0开始的而matlab中的数组的下标是从1开始的，又比如说C+里面的变量都必须先定义之后才能使用，函数的调用以及要用到指向函数的指针，在语法规则的允许范围内，对函数以及指针

13、也可以灵活的使用，并且我们应该要理解函数调用的实质，即实参与形参所共用同一段代码，即将实参的地址传递给形参，而C+程序里面的函数名，数组名等都可以代表函数以及数组的首地址，因此，函数调用时只需要将对应的数组名或者变量名放入调用函数中即可以实现函数的调用。当然，我所学到的是远远不够的，不仅仅是学习方面，在其他方面也应该努力去学习，课程设计只是给了我们一个训练和发掘自己的一个平台，最为重要的是在于我们以怎么的一种态度去学习，而不仅仅是把他当成一种必须的任务去完成，我们对于新事物和不懂得东西应该充满着激情和好奇心，不断培养自己的兴趣和爱好。学习东西要有信心和决心，让生活的每一天都过得很充实，这应该成

14、为每一个人的梦想。参 考 文 献 1谭浩强编.C+ 程序设计 M.北京：清华大学出版社，20052谭浩强主编.C程序设计题解与上机指导（第二版） M.北京：清华大学出版社，20053John H.Mathews Kurtis D.Fink 著 数值方法（MATLAB版）（第四版）M.北京：电子工业出版社，2002附录1用连续梯形公式求解积分程序代码#include#includeusing namespace std;/被积分函数定义double f(double x)return(x*x*x);/定义结束double tt(double(*ff)(double x),double a,dou

15、ble b,double eps)double T100,H100,s=0; int i,n=0; T0=(b-a)*(ff(a)+ff(b)/2;H0=(b-a)*(ff(0.5*(a+b);docoutTn Hnendl;/ 每次循环都输出Tn，Hn的值Tn+1=0.5*(Tn+Hn);s=0; for(i=1;i=eps);/循环执行条件判断coutTn Hna1b1eps;ff=f;result=tt(ff,a1,b1,eps);/函数调用coutresult=resultendl;return 0;2高斯选主元解线性方程组程序代码#include #define m 10 #incl

16、udeusing namespace std;double Amm,bm,Augmm+1,s;/定义全局变量void equ(double Amm,double bm,double xm,int n) int i,i1,j,k; double Augmm+1,max,temp,l;/构造Aug for(i=0;in;i+) /找主元 for(j=0;jn;j+) Augij=Aij; Augin=bi; for(i1=0;i1n-1;i1+) max=fabs(Augi1i1);k=i1; for(i=i1;in;i+) if(maxfabs(Augii1) max=fabs(Augii1);

17、 k=i; /换行 for(j=i1;jn+1;j+) temp=Augi1j; Augi1j=Augkj; Augkj=temp; /高斯消元，以第i1行为工具行处理一下各行 for (k=i1+1;kn;k+) l=-Augki1/Augi1i1; for(j=i1;j=0;i-) s=0; for(j=i+1;jn;j+) s=s+Augij*xj; xi=(Augin-s)/Augii; int main() int i,j,n; double Amm,bm,xm; cout请输入未知量个数nn; if(nm) cout问题规模太大，请重新定义符号常量m的值endl; for(i=0;

18、in;i+) cout请输入A的i+1行:; for(j=0;jAij; cout请输入b:; for(j=0;jbj; equ(A,b,x,n);cout输出结果:x=;for(i=0;in;i+) coutxi ;cout;endl;return 0; 3经典四阶龙格库塔法解一阶微分方程组程序代码#include #include #define N 20using namespace std;/定义微分方程float fx(float dd,float pp) float der; der=dd+pp;return(der);float fy(float d,float p) float

19、 derivat; derivat=-d-p;return(derivat);void main()float t0,t,h,nn; int n; /用于对nn取整 float k1,k2,k3,k4,k11,k22,k33,k44; float xN,yN; /存放计算出的微分方程的数值解float xx,yy; /精确值 float errorx,errory;/误差int i=0,j;coutt0;coutt;coutx0;couty0;couth; nn=(t-t0)/h; coutnn=nnendl;n=int(nn); /强制类型转换 coutn=nendl; for(j=0;j=

20、n;j+) xx=3*exp(t0-1)-t0-1; yy=exp(1-t0)-t0+1; errorx=xi-xx; errory=yi-yy; coutxt0=xitxx=xxterrorx=errorx; coutxt0=yityy=yyterrory=errory; k1=h*fx(t0,xi); /求K1 k2=h*fx(t0+h/2),(xi+k1*h/2); /求K2 k3=h*fx(t0+h/2),(xi+k2*h/2); /求K3 k4=h*fx(t0+h),(xi+k3*h); /求K4 xi+1=xi+(k1+2*k2+2*k3+k4)/6); /求xi+1 k11=h*

21、fy(t0,yi); /求K11 k22=h*fy(t0+h/2),(yi+k11*h/2); /求K22 k33=h*fy(t0+h/2),(yi+k22*h/2); /求K33 k44=h*fy(t0+h),(yi+k33*h); /求K44 yi+1=yi+(k11+2*k22+2*k33+k44)/6); /求yi+1 t0=t0+float(h); i+;4三次样条插值算法（紧压样条）程序代码#includeusing namespace std;#define MAX_N 20typedef struct tagPOINT/定义一个结构体类型 double x;double y;P

22、OINT;int main() int n,i,k;POINT pointsMAX_N+1;double hMAX_N+1,bMAX_N+1,cMAX_N+1,dMAX_N+1,MMAX_N+1;double uMAX_N+1,vMAX_N+1,yMAX_N+1;double x,p,q,S;coutinput n value:n;if (nMAX_N)cout输入过大，请重新输入endl;return 1;if(n=0) cout 请重新输入MAX_Nendl;return 1;coutNow input the (x_i,y_i),i=0.nendl;for(i=0;ipointsi.xp

23、ointsi.y;coutM0;coutMn;coutx;if(xpointsn.x|xpoints0.x) cout请输入一个数从points0.x到pointsn.xendl;return 1;h0=points1.x-points0.x;for(i=1;in;i+)hi=pointsi+1.x-pointsi.x;bi=hi/(hi+hi-1);ci=1-bi;di=6*(pointsi+1.y-pointsi.y)/hi-(pointsi.y-pointsi-1.y)/hi-1)/(hi+hi-1);d1-=c1*M0;dn-1-=bn-1*Mn;bn-1=0;c1=0;v0=0;fo

24、r (i=1;in;i+)ui=2-ci*vi-1;vi=bi/ui;yi=(di-ci*yi-1)/ui;for(i=1;i=pointsk.x)k+;k=k-1;p=pointsk+1.x-x;q=x-pointsk.x;S=(p*p*p*Mk+q*q*q*Mk+1)/(6*hk)+(p*pointsk.y+q*pointsk+1.y)/hk-hk*(p*Mk+q*Mk+1)/6;coutS(x)=Sendl;5龙贝格求积分算法程序代码#includeusing namespace std;#include#define f(x) (cos(x)#define N_H 20#define

25、MAXREPT 10#define eps 0.00001#define a 1#define b 2/定义函数double computeT(double aa,double bb,int n)int i;/将非端点处的节点处的函数值不断累积求和再加上端点函数值和的平均值，/最后返回和与步长的乘积double sum,h=(bb-aa)/n;sum=0;/在每次循环之前讲sum清零for(i=1;in;i+)sum=sum+f(aa+i*h);sum=sum+(f(aa)+f(bb)/2;return(h*sum);void main()int i;int n=N_H,m=0;double

26、T2MAXREPT+1;/定义二维数组T10=computeT(a,b,n);n=2*n;for(m=1;mMAXREPT;m+)for(i=0;im;i+)T0i=T1i;T10=computeT(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)eps)coutT=T1mendl;return;也可以直接定义被积函数，然后利用函数指针实现调用#includeusing namespace std;#include#define N_H 20#define MAXREPT 10#d

27、efine eps 0.00001float a=1,b=2; /a,b为积分区间/定义被积函数double f(double x)return(cos(x);double computeT(double(*ff)(double x),double aa,double bb,int n)int i;double sum,h=(bb-aa)/n;ff=f;sum=0; /在每次循环之前讲sum清零for(i=1;in;i+)sum=sum+f(aa+i*h);sum=sum+(f(aa)+f(bb)/2;return(h*sum);void main()int i;int n=N_H,m=0;d

28、ouble T2MAXREPT+1;/定义二维数组T10=computeT(f,a,b,n);n=2*n;for(m=1;mMAXREPT;m+)for(i=0;im;i+)T0i=T1i;T10=computeT(f,a,b,n);/函数调用，此时的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)eps)coutT=T1mendl;return;6 M次多项式曲线拟合程序代码#include#include#define m 10#define n 50using name

29、space std;double Amm,bm,Augmm+1,s;/定义全局变量void equ(double Am+1m+1,double bm+1,double xm+1,int n2) int i,i1,j,k; double Augm+1m+2,max,temp,l;/构造Aug for(i=0;in2;i+) for(j=0;jn2;j+) Augij=Aij; Augin2=bi; for(i1=0;i1n2-1;i1+) max=fabs(Augi1i1);k=i1; for(i=i1;in2;i+) if(maxfabs(Augii1) max=fabs(Augii1); k

30、=i; /换行 for(j=i1;jn2+1;j+) temp=Augi1j; Augi1j=Augkj; Augkj=temp; /高斯消元，以第i1行为工具行处理一下各行 for (k=i1+1;kn2;k+) l=-Augki1/Augi1i1; for(j=i1;j=0;i-) s=0;/每次循环都将s清零，以防后一次循环时出错 for(j=i+1;jn2;j+) s=s+Augij*xj; xi=(Augin2-s)/Augii; void polfit(double xn,double yn,double am+1,int m1,int n1)double Am+1m+1,bm+1

31、,s;int i,k,p;for(k=0;km1+1;k+)for(p=0;pm1+1;p+)/矩阵中各个幂的和s=0;for(i=0;in1;i+)s=s+pow(xi,k+p);Akp=s; /A的m+1列放入常数列向量s=0;for(i=0;in1;i+)s=s+pow(xi,k)*yi;bk=s;/求解Aa=bequ(A,b,a,m1+1);int main()int i,m1,n1;double xn,yn,a1m+1;coutn1;coutendl;coutm1;cout输入x的元素:;for(i=0;ixi;coutendl;cout输入y的元素:;for(i=0;iyi;pol

32、fit(x,y,a1,m1,n1);cout拟合多项式的系数为:(按照升幂排列);for(i=0;im1+1;i+)couta1i ;coutendl;couty=a10+(a11x)+a12x2endl;return 0;7牛顿 拉夫森迭代法程序代码#includeusing namespace std;#include/定义函数#define f(x) (x*x-3)*x-1)/通过xk+1=xk+f(xk)/f(xk)构造迭代式iterate(x)#define iterate(x) (x-(x*x-3)*x-1)/(3*x*x-3)void main()double eps=0.000

33、1,x0=1.5;/给定初始迭代点及精度要求int i;double xk=x0,xk1=x0;for(i=0;i1000;i+) cout第i+1次迭代得：xk1endl;xk1=iterate(xk);if (fabs(xk1-xk)eps|fabs(f(xk1)eps) cout!Root :xk1endl;return;xk=xk1;8牛顿法解非线性方程组程序代码#include#include#include#define N 2 #define epsilon 0.0001 #define Max 100 using namespace std;const int N2=2*N;void ff(float xxN,float yyN)float x,y; int i; x=xx0;y=xx1;yy0=x*x-y-0.2;yy1=y*y-x-0.3; cout当前自变量向量为：endl;for( i=0;iN;i+) coutsetw(10)setiosflags(ios:right)setprecision(5)setiosflags(ios:fixed)xxi ; coutendl;coutendl; cout当前因变量向量为：endl;for( i=0;iN;i+)