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

数值计算课程设计典型数值算法的C+语言程序设计1.经典四阶龙格库塔法解一阶微分方程组1.1算法说明龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。4阶龙格-库塔方法(RK4)可模拟N=4的泰勒方法的精度。这种算法可以描述为，自初始点开始，利用下面的计算方法生成近似序列1.2经典四阶龙格库塔法解一阶微分方程组算法流程图图1-1算法流程图1.3经典四阶龙格库塔法解一阶微分方程组程序调试将编写好的代码放在VC6.0环境中编译，直接执行程序便可以得到求解微分方程，并且的结果。如图：图1-2 程序调试图将这些点进行插值或者拟合后就可以得到微分方程的解。1.4 经典四阶龙格库塔法解一阶微分方程组代码#include #include using namespace std;/f为函数的入口地址，x0、y0为初值，xn为所求点，step为计算次数double Runge_Kuta( double (*f)(double x, double y), double x0, double y0, double xn, long step ) double k1,k2,k3,k4,result; double h=(xn-x0)/step; if(step=0) return(y0); if(step=1) k1=f(x0,y0); k2=f(x0+h/2, y0+h*k1/2); k3=f(x0+h/2, y0+h*k2/2); k4=f(x0+h, y0+h*k3); result=y0+h*(k1+2*k2+2*k3+k4)/6; else double x1,y1; x1=xn-h; y1=Runge_Kuta(f, x0, y0, xn-h,step-1); k1=f(x1,y1); k2=f(x1+h/2, y1+h*k1/2); k3=f(x1+h/2, y1+h*k2/2); k4=f(x1+h, y1+h*k3); result=y1+h*(k1+2*k2+2*k3+k4)/6; return(result);int main() double f(double x, double y); double x0,y0; double a,b; coutx0y0; coutab; /double x0=0,y0=1; double x,y,step; long i; coutstep; /step=0.1; cout.precision(10); for(i=0;i=(b-a)/step;i+) x=x0+i*step; coutsetw(8)xsetw(18)Runge_Kuta(f,x0,y0,x,i)endl; return(0);double f(double x, double y)35数值计算课程设计 double r; r=(x-y)/2; return(r);2. 高斯列主元法解线性方程组2.1算法说明首先将线性方程组做成增光矩阵，对增广矩阵进行行变换。对第元素，在第i列中，第i行及以下的元素选取绝对值最大的元素，将该元素最大的行与第i行交换，然后采用高斯消元法将新得到的消去第i行以下的元素。一次进行直到。从而得到上三角矩阵。再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。2.2高斯列主元算法流程图图2-1 算法流程图2.3高斯列主元程序调试对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数为4，输入4行5列的增广矩阵：图2-1 程序调试图1按回车键后，程序执行得如下所示的结果：图2-2 程序调试图22.4 高斯列主元算法代码#include#include#include using namespace std;/在列向量中寻找绝对值最大的项，并返回该项的标号int FindMax(int p,int N,double *A)int i=0,j=0;double max=0.0;for(i=p;imax)j=i;max=fabs(Ai*(N+1)+p);return j;/交换矩阵中的两行void ExchangeRow(int p,int j,double *A,int N)int i=0;double C=0.0;for(i=0;iN+1;i+)C=Ap*(N+1)+i;Ap*(N+1)+i=Aj*(N+1)+i;Aj*(N+1)+i=C;/上三角变换，A为增广矩阵的指针，N为矩阵的行数。void uptrbk(double *A,int N)int p=0,k=0,q=0,j=0;double m=0.0;for(p=0;pN-1;p+)/找出该列最大项的标号j=FindMax(p,N,A);/交换p行和j行ExchangeRow(p,j,A,N);if(Ap*(N+1)+p=0)printf(矩阵是一个奇异矩阵。没有唯一解。);break;/消去P元素一下的p列内容。for(k=p+1;kN;k+)m=Ak*(N+1)+p/Ap*(N+1)+p;for(q=p;qN+1;q+)Ak*(N+1)+q=Ak*(N+1)+q-m*Ap*(N+1)+q; coutn增广矩阵高斯列主元消去后的矩阵为:n;for(j=0;jN*(N+1);j+)if(j%(N+1)=0)coutendl;cout=0;k-)temp=0.0;for(i=k+1;iN;i+)temp=temp+Ak*(N+1)+i*Xi;Xk=(Ak*(N+1)+N-temp)/Ak*(N+1)+k;return X;main()int N=0,i=0;double *A=NULL,*X=NULL; coutN;A=(double*)calloc(N*(N+1),sizeof(double); cout请输入待求解方程组的增广矩阵(%d行%d列):n,N,N+1;for(i=0;iAi;system(cls);cout方程的增广矩阵为:n;for(i=0;iN*(N+1);i+)if(i%(N+1)=0)printf(n);printf(%lft,Ai);uptrbk(A,N);X=backsub(A,N);coutnn方程组的解为:n;for(i=0;iN;i+)coutX()=Xi;free(A);free(X);exit(0);3.牛顿法解非线性方程组3.1算法说明设已知。第1步：计算函数第2步：计算雅可比矩阵第3步：求线性方程组的解。第4步：计算下一点重复上述过程。3.2 牛顿法解非线性方程组算法流程图图3-1 算法流程图3.3 牛顿法解非线性方程组算法程序调试我们以方程组为例进行求解，我们按命令的要求，依次输入牛顿法解非线性方程组:x2-2*x-y+0.5=0x2+4*y2-4=0输入的初始近似值x0,y02.00 0.25请依次输入P的误差限，F(P)的误差限，最大迭代次数0.0001 0.0001 1000收敛到P的解为:X(1)=1.900691X(2)=0.311213迭代次数为:2误差为:0.0000003.4牛顿法解非线性方程组算法程序代码#include#include#includeusing namespace std;#define N 2 /用来设置方程组的行数#define eps 2.2204e-16double* MatrixMultiply(double* J,double Y);double* Inv(double *J);double norm(double Q);double* F(double X);double* JF(double X);int method(double* Y,double epsilon);int newdim(double P,double delta,double epsilon,int max1,double *err)double *Y=NULL,*J=NULL,Q2=0,*Z=NULL,*temp=NULL;double relerr=0.0;int k=0,i=0,iter=0;Y=F(P);for(k=1;kmax1;k+)J=JF(P);temp=MatrixMultiply(Inv(J),Y);for(i=0;iN;i+)Qi=Pi-tempi;Z=F(Q);for(i=0;iN;i+)tempi=Qi-Pi;*err=norm(temp);relerr=*err/(norm(Q)+eps);for(i=0;iN;i+)Pi=Qi;for(i=0;iN;i+)Yi=Zi;iter=k;if(*errdelta)|(relerrdelta)|method(Y,epsilon)break;return iter;int method(double* Y,double epsilon)if(fabs(Y0)epsilon&fabs(Y1)epsilon)return 1;elsereturn 0;/矩阵乘法，要求，J为方阵，Y为与J维数相同的列向量double *MatrixMultiply(double* J,double Y)double *X=NULL;int i=0,j=0;X=(double*)malloc(N*sizeof(double);for(i=0;iN;i+)Xi=0;for(i=0;iN;i+)for(j=0;jN;j+)Xi+=Ji*N+j*Yj;return X;/二阶矩阵的求逆（在M次多项式曲线拟合算法文件中给出了对任意可逆矩阵的求逆算法）double *Inv(double *J)double X4=0,temp=0.0;int i=0;temp=1/(J0*J3-J1*J2);X0=J3;X1=-J1;X2=-J2;X3=J0;for(i=0;i4;i+)Ji=temp*Xi;return J;double norm(double Q)double max=0.0;int i=0;for(i=0;imax)max=Qi;return max;double* F(double X)double x=X0;double y=X1;double *Z=NULL;Z=(double*)malloc(2*sizeof(double);Z0=x*x-2*x-y+0.5;Z1=x*x+4*y*y-4;return Z;double* JF(double X)double x=X0;double y=X1;double *W=NULL;W=(double*)malloc(4*sizeof(double);W0=2*x-2;W1=-1;W2=2*x;W3=8*y;return W; main()double P2=0;double delta=0.0,epsilon=0.0,err=0.0; int max1=0,iter=0,i=0;cout牛顿法解非线性方程组:nx2-2*x-y+0.5=0nx2+4*y2-4=0n; coutn输入的初始近似值x0,y0n;for(i=0;iPi;coutdelta;cinepsilon;cinerr;iter=newdim(P,delta,epsilon,max1,&err);cout收敛到P的解为:n;for(i=0;i2;i+)coutX(i+1)=Piendl;coutn迭代次数为: iter;coutn误差为:errendl;getchar();4.龙贝格求积分算法4.1算法说明生成的逼近表，并以为最终解来逼近积分逼近存在于一个特别的下三角矩阵中，第0列元素用基于个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为当时，程序在第行结束。4.1龙贝格求积分算法流程图图4-1 算法流程图4.3龙贝格求积分算法程序调试我们以求解积分方程为例，对所编写的龙贝格求积分算法程序进行编译和链接，经执行后得如下所示的窗口图4-2 程序调试图4.4龙贝格求积分算法代码#include#includeusing namespace std;double f(double x)returnx*x;main()int M=1,n=0,p=0,K=0,i=0,j=0,J=0;double h=0.0,a=0.0,b=0.0,err=1.0,quad=0.0,s=0.0,x=0.0,tol=0.0;double R3030=0;a=0;b=1;h=b-a;n=4;tol=0.000001;cout求解函数y=x*x在(0,1)上的龙贝格矩阵endl;cout龙贝格矩阵最大行数为nendl;cout误差限为toltol)&(Jn)|(J4)J=J+1;h=h/2;s=0;for(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(K=1;K=J;K+)RJK=RJK-1+(RJK-1-RJ-1K-1)/(pow(4,K)-1);err=fabs(RJ-1J-1-RJK);quad=RJJ;coutn龙贝格矩阵为:n;for(i=0;i(J+1);i+)for(j=0;j(J+1);j+)printf(%.5lf ,Rij);printf(n);coutn积分值为:quad=quad;coutn误差估计为:err=%lf,err;coutn使用过的最小步长:h=hendl;getchar();5.三次样条插值算法5.1算法说明策略描述包含和的方程(i)三次紧压样条，确定，（如果导数已知，这是“最佳选择”）(ii)natural三次样条（一条“松弛曲线”），(iii)外挂到端点(iv) 是靠近端点的常量，(v)在每个端点处指定，5.2 三次样条插值算法（压紧样条）程序调试求三次紧压样条曲线，我们以经过点（0,0），（1,0.5）,(2,2.0)和（3,1.5），且一阶导数的边界条件为S(0)=0.2和S(3)=-1为例。我们将所编写的程序经过编译，链接和执行后得如下所示的结果图5-1 程序调试图我们借助Matlab绘制出以上三次压紧样条的函数图像如下所示图5-2 三次压紧样条的函数图像5.3 三次样条插值算法（压紧样条）代码#include#includeusing namespace std;const MAX= 4;double *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;cout求解经过点（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-1;M0=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;cout求得的三次压紧样条曲线的矩阵S为:n;for(i=0;iMAX-1;i+)for(k=0;kMAX;k+)couttSik;coutendl;for(i=0;iN;i+) coutn在区间（0,1)上的样条为：y=+Si0x3+Si1x2+Si2x+Si3;/printf(n在区间(0,1)上的样条为:y=%+lfx3%+lfx2%+lfx%+lf,Si0,Si1,Si2,Si3);getchar();6.M次多项式曲线拟合6.1算法说明设有N个点，横坐标是确定的。最小二乘抛物线的系数表示为求解A，B和C的线性方程组为6.2 M次多项式曲线拟合算法流程图图6-1 算法流程图6.3 M次多项式曲线拟合算法程序调试我们按命令依次输入命令如下命令后，得程序执行结果如下图6-2 程序调试图6.4 M次多项式曲线拟合算法代码#include#include#include using namespace std;const MAX =20;/求解任意可逆矩阵的逆，X为待求解矩阵，E为全零矩阵，非单位矩阵，也可以是单位矩阵void inv(double XMAXMAX,int n,double EMAXMAX)int i=0,j=0,k=0;double temp=0.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;main()int n=0,M=0,i=0,j=0,k=0;double XMAX=0,YMAX=0,FMAXMAX=0,BMAX=0;double AMAXMAX=0,BFMAXMAX=0,EMAXMAX=0,CMAX=0;coutn;coutn请输入n个点的X坐标序列:n;for(i=0;iXi;coutn请输入n个点的Y坐标序列:n;for(i=0;iYi;coutM;for(i=0;in;i+)for(k=1;k=M+1;k+)Fik-1=pow(Xi,k-1);/求F的转置for(i=0;in;i+)for(j=0;jM+1;j+)BFji=Fij; /计算F与其转置的BF的乘 for(i=0;iM+1;i+)for(j=0;jM+1;j+)for(k=0;kn;k+)Aij+=BFik*Fkj;/计算F的转置BF与Y的乘for(i=0;iM+1;i+)for(j=0;jn;j+)Bi+=BFij*Yj;inv(A,n,E);/计算A的逆E与B的乘for(i=0;iM+1;i+)for(j=0;jn;j+)Ci+=Eij*Bj;coutn拟合后的M=0;i-)coutCit;coutn拟合后的M次多项式为:n;cout=0;i-)if(i=0)cout+setprecision(3)Ciendl;/printf(%+.3lfn,Ci);elsecout+setprecision(3)Ci*xiendl;getchar();7.二分法解非线性方程7.1算法说明1.是起始区间，是中点。2.是第二个区间，它包含零点r，同时是中点，区间的宽度范围是的一半。3.得到第n个区间（包含r，并有中点）后，可构造出，它也包括r，宽度范围是的一半。7.2 二分法解非线性方程算法流程图图7-1 算法流程图7.3 二分法解非线性方程算法程序调试 我们将所编写的二分法算法程序经过编译，链接和执行后得所下结果图7-2 程序调试图7.4 二分法解非线性方程算法代码#includeusing namespace std;const eps =1e-10;double f(double x)return 3*x*x*x+5; void main()double ga=0.0,gb=0.0,gc=0.0,a=0.0,b=0.0,c=0.0;cout用二分法寻找方程3*x3+5=0的根endl;cout求根区间为（-5，5）eps)c=(a+b)/2;gc=f(c);if(gc=0)break;else if(gc*gb0)a=c;ga=gc;elseb=c;gb=gc;cout方程的根为:X=bendl;getchar();8.不动点法解非线性方程8.1算法说明先将改写成然后对进行迭代，即 其中然后判断是否成立，成立则返回，不成立就重复以上步骤8.2 不动点法解非线性方程算法程序调试我们将编写好的不咪法解非线性方程算法程序进行编译，链接和执行后得如下所示结果图8-1 程序调试图8.3 不动点法解非线性方程算法代码#include#includeusing namespace std;const MAX= 20;const eps =2.2204e-16;double g(double x)return 1-x*x/4+x;main()double PMAX=0,err=0.0,relerr=0.0,tol=0.0,p=0.0,p0=0.0;int k=0,max1=0,i=0;cout不动点法解非线性方程f(x)=1-x2/2n;cout方程在0，1上有解，初始值为p0=0n;/初始化P0=p0=0;max1=100;tol=0.001;for(k=2;k=max1;k+) Pk-1=g(Pk-2);err=fabs(Pk-1-Pk-2);relerr=err/(fabs(Pk-1+eps);p=Pk-1;if(errtol)|(relerrtol)break;if(k=max1)cout迭代次数超过允许的最大迭代次数！;coutn不动点的近似值为:p; coutn程序迭代次数为:k;coutn近似值的误差为:err;printf(n求解不动点近似值的序列:n);for(i=0;ik;i+)coutPi;getchar();9.霍纳法多项式求值9.1算法说明设是按式(20)给出的多项式，且是用于计算的数。设，并计算 其中则。进一步考虑，如果则其中是n-1阶多项式的商，是余数。9.2霍纳法多项式求值算法流程图图9-1算法流程图9.3霍纳法多项式求值算法程序调试我们将所编写的霍纳法多项式求值算法程序经过编译，链接和执行，根据窗口的命令依次输入如下的命令后，得到如下所示的结果图9-2 程序调试图9.4霍纳法多项式求值算法代码#include#include#includeusing namespace std;void main()int i=0,n=0;double *A=NULL,*B=NULL,c=0.0;cout霍纳法求解多项式endl;coutn;A=(double*)malloc(n*sizeof(double);B=(double*)malloc(n*sizeof(double);cout=0;i-)cinAi;coutc;Bn-1=An-1;for(i=n-2;i=0;i-)Bi=Ai+c*Bi+1;cout=0;i-)if(i=0)cout+setprecision(2)Aiendl;elsecout+setprecision(2)Ai*xiendl;coutn当x=c时的解为P(x)=B0endl; free(A);free(B);getchar();10. 牛顿-拉弗森迭代解非线性方程10.1算法说明使用初始近似值，利用迭代 其中计算函数的根的近似值。10.2牛顿-拉弗森迭代解非线性方程算法程序调试 我们将所编写的牛顿-拉弗森迭代解非线性方程算法程序经过编译，链接和执行后得如下所示的结果图10-1程序调试图10.3牛顿-拉弗森迭代解非线性方程算法代码#include#includeusing namespace std;double f(double x)return x*x-2*x-1;double df(double x)return 2*x-2;void mai