计算方法程序(C++)

1、1.经典四阶龙格库塔法解一阶微分方程组代码#include<stdio.h>#include<stdlib.h>#define N 24/4阶龙格-库塔方法void rk4(double (*p)(double t,double y),double a,double b,double ya,int M)double h=0.0,k1=0.0,k2=0.0,k3=0.0,k4=0.0;double *T=NULL,*Y=NULL;double *R2;int i=0,j=0;h=(b-a)/M;printf("h=%fnn",h);T=(double*

2、)malloc(M+1)*sizeof(double);/T=0;Y=(double*)malloc(M+1)*sizeof(double);/T=0;T0=a;for(i=1;i<M+1;i+)Ti=Ti-1+h;Y0=ya;for(i=0;i<M;i+)k1=h*p(Ti,Yi);k2=h*p(Ti+h/2,Yi+k1/2);k3=h*p(Ti+h/2,Yi+k2/2);k4=h*p(Ti+h,Yi+k3);Yi+1=Yi+(k1+2*k2+2*k3+k4)/6;R0=T;R1=Y;printf("tk=ty(tk)=n");for(i=0;i<N+1

3、;i+)for(j=0;j<2;j+)printf("%lf ",Rji);printf("n");/return R;/微分方程的函数文件double f(double t,double y)return (t-y)/2;main()int i=0,j=0;printf("用4阶龙格-库塔方法解微分方程：y'=(t-y)/2, y(0)=1n");printf("微分方程的解为:n");rk4(f,0,3,1.0,N);return 0;2.高斯列主元算法代码#include<stdio.h&

4、gt;#include<stdlib.h>#include<math.h>/在列向量中寻找绝对值最大的项，并返回该项的标号int FindMax(int p,int N,double *A)int i=0,j=0;double max=0.0;for(i=p;i<N;i+)if(fabs(Ai*(N+1)+p)>max)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;i<

5、;N+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;p<N-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+

6、1;k<N;k+)m=Ak*(N+1)+p/Ap*(N+1)+p;for(q=p;q<N+1;q+)Ak*(N+1)+q=Ak*(N+1)+q-m*Ap*(N+1)+q;printf("n增广矩阵高斯列主元消去后的矩阵为:n");for(j=0;j<N*(N+1);j+)if(j%(N+1)=0)printf("n");printf("%lft",Aj);/下面是回代函数double* backsub(double *A,int N)double* X=NULL,temp=0.0;int k=0,i=0;X=(dou

7、ble*)malloc(N*sizeof(double);XN-1=A(N-1)*(N+1)+N/A(N-1)*(N+1)+N-1;for(k=N-2;k>=0;k-)temp=0.0;for(i=k+1;i<N;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;printf("n请输入待求解方程组的增广矩阵的行数:");scanf("%d",&N);A=(double

8、*)calloc(N*(N+1),sizeof(double);printf("请输入待求解方程组的增广矩阵(%d行%d列):n",N,N+1);for(i=0;i<N*(N+1);i+)scanf("%lf",&Ai);system("cls");printf("方程的增广矩阵为:n");for(i=0;i<N*(N+1);i+)if(i%(N+1)=0)printf("n");printf("%lft",Ai);uptrbk(A,N);X=backsu

9、b(A,N);printf("nn方程组的解为:n");for(i=0;i<N;i+)printf("X(%d)= %lfn",i+1,Xi);free(A);free(X);exit(0);3.牛顿法解非线性方程组算法程序代码#include<stdio.h>#include<stdlib.h>#include<math.h>#define N 2 /用来设置方程组的行数#define eps 2.2204e-16double* MatrixMultiply(double* J,double Y);double

10、* 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;k<max1;k+)J=

11、JF(P);temp=MatrixMultiply(Inv(J),Y);for(i=0;i<N;i+)Qi=Pi-tempi;Z=F(Q);for(i=0;i<N;i+)tempi=Qi-Pi;*err=norm(temp);relerr=*err/(norm(Q)+eps);for(i=0;i<N;i+)Pi=Qi;for(i=0;i<N;i+)Yi=Zi;iter=k;if(*err<delta)|(relerr<delta)|method(Y,epsilon)break;return iter;int method(double* Y,double e

12、psilon)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;i<N;i+)Xi=0;for(i=0;i<N;i+)for(j=0;j<N;j+)Xi+=Ji*N+j*Yj;return X;/二阶矩阵的求

13、逆（在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;i<4;i+)Ji=temp*Xi;return J;double norm(double Q)double max=0.0;int i=0;for(i=0;i<N;i+)if(Qi>max)max=Qi;return max;double* F(double X)double x=X0;double

14、y=X1;double *Z=NULL;Z=(double*)malloc(2*sizeof(double);Z0=x*x-2*x-4*y+4;Z1=x*x+2*y*y-5;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=

15、0,i=0;printf("牛顿法解非线性方程组:nx*x-2*x-4*y+4=0nx*x+2*y*y-5=0n");printf("n输入的初始近似值(建议使用2.00 0.25)n");for(i=0;i<2;i+)scanf("%lf",&Pi);printf("请依次输入P的误差限，F(P)的误差限，最大迭代次数n");scanf("%lf%lf%d",&delta,&epsilon,&max1);iter=newdim(P,delta,epsilo

16、n,max1,&err);printf("收敛到P的解为:n");for(i=0;i<2;i+)printf("X(%d)=%lfn",i+1,Pi);printf("n迭代次数为:%d",iter);printf("n误差为:%lfn",err);return 0;4.龙贝格求积分算法代码#include<stdio.h>#include<math.h>double f(double x)returnx*x;main()int M=1,n=0,p=0,K=0,i=0,j=0,

17、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;printf("求解函数y=x*x在(0,1)上的龙贝格矩阵n");printf("龙贝格矩阵最大行数为%dn误差限为%lfn",n,tol);R00=h*(f(a)+f(b)/2;while(err>tol)&&(J<n)|(J<4)J=J+1;h=h/2;s=0;for(p=1;p<=M;p

18、+)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;printf("n龙贝格矩阵为:n");for(i=0;i<(J+1);i+)for(j=0;j<(J+1);j+)printf("%.5lf ",Rij);printf("n");printf("n积分值为:quad=%lf",qua

19、d);printf("n误差估计为:err=%lf",err);printf("n使用过的最小步长:h=%lfn",h);getch();5. 三次样条插值算法（压紧样条）代码#include<stdio.h>#include<stdlib.h>#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;

20、double *divide(double Y,int N,double H)int i=0;double *D=NULL;D=(double*)malloc(N*sizeof(double);for(i=0;i<N;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=NUL

21、L,*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;i<N-2;i+)Ai=Hi+1;for(i=0;i<N-1;i+)Bi=2*(Hi+Hi+1);for(i=1;i<N-1;i+)Ci=Hi+1;U=diff(D,N);for(i=0;i<N;i+)Ui=Ui*6;B0=B0-H0

22、/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<=N-1;k+) temp=Ak-2/Bk-2; Bk-1=Bk-1-temp*Ck-2; Uk-1=Uk-1-temp*Uk-2;MN-1=UN-2/BN-2;for(k=N-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

23、/2;Sk2=Dk-Hk*(2*Mk+Mk+1)/6;Sk3=Yk;printf("求得的三次压紧样条曲线的矩阵S为:n");for(i=0;i<MAX-1;i+)for(k=0;k<MAX;k+)printf("%lft",Sik);printf("n");for(i=0;i<N;i+)printf("n在区间(0,1)上的样条为:y=%+lfx3%+lfx2%+lfx%+lf",Si0,Si1,Si2,Si3);return 0;6. M次多项式曲线拟合算法代码#include<stdi

24、o.h>#include<math.h>#define 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;i<MAX;i+)for(j=0;j<MAX;j+)if(i=j)Eij=1;for(i=0;i<n-1;i+)temp=Xii;for(j=0;j<n;j+)Xij=Xij/temp;Eij=Eij/temp;for(k=0;k<

25、n;k+)if(k=i)continue;temp=-Xii*Xki;for(j=0;j<n;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;printf("tttM次多项式曲线拟合nn请先输入待拟合的点的个数:");scanf("%d",&n);printf("n请输入%d个点的X

26、坐标序列:n",n);for(i=0;i<n;i+)scanf("%lf",&Xi);printf("n请输入%d个点的Y坐标序列:n",n);for(i=0;i<n;i+)scanf("%lf",&Yi);printf("n请输入需要拟合的次数:");scanf("%d",&M);for(i=0;i<n;i+)for(k=1;k<=M+1;k+)Fik-1=pow(Xi,k-1);/求F的转置for(i=0;i<n;i+)for(

27、j=0;j<M+1;j+)BFji=Fij; /计算F与其转置的BF的乘for(i=0;i<M+1;i+)for(j=0;j<M+1;j+)for(k=0;k<n;k+)Aij+=BFik*Fkj;/计算F的转置BF与Y的乘for(i=0;i<M+1;i+)for(j=0;j<n;j+)Bi+=BFij*Yj;inv(A,n,E);/计算A的逆E与B的乘for(i=0;i<M+1;i+)for(j=0;j<n;j+)Ci+=Eij*Bj;printf("n拟合后的%d次多项式系数为，幂次由高到低：n",M);for(i=M;i

28、>=0;i-)printf("%lft",Ci);printf("n拟合后的%d次多项式为:n",M);printf("nP(x)=");for(i=M;i>=0;i-)if(i=0)printf("%+.3lfn",Ci);elseprintf("%+.3lf*x%d",Ci,i);return 0;7. 不动点法解非线性方程算法代码#include<stdio.h>#include<math.h>#define MAX 20#define eps 2.22

29、04e-16double g(double x)return 1-x*x/9+2*x+3;void printx()int i=0;for(i=0;i<20;i+)printf(" ");for(i=0;i<40;i+)printf("*");for(i=0;i<20;i+)printf(" ");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;printf("nnnn");printx(

30、);printf("ttt不动点法解非线性方程f(x)= 1-x*x/9+2*x+3n");printf("ttt方程在0，10上有解，初始值为p0=0n");printx();/初始化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(err<tol)|(relerr<tol)break;if(k=max1)printf("迭代次数超过允许的最大

31、迭代次数！");printf("nnnttt不动点的近似值为:%lf",p);printf("nttt程序迭代次数为%d",k);printf("nttt近似值的误差为:%lf",err);printf("nttt求解不动点近似值的序列:n");for(i=0;i<k;i+)printf("%lf ",Pi);return 0;8.二分法解非线性方程算法代码#include<stdio.h>#define eps 1e-10double f(double x)retu

32、rn 6*x*x*x+19*x+38;void printx()int i=0;for(i=0;i<20;i+)printf(" ");for(i=0;i<40;i+)printf("*");for(i=0;i<20;i+)printf(" ");main()double ga=0.0,gb=0.0,gc=0.0,a=0.0,b=0.0,c=0.0;printf("nnnnn");printx();printf("ttt用二分法寻找方程6*x3+19*x+38=0的根n");p

33、rintf("ttt求根区间为（-5，5）n");a=-5,b=5;ga=f(a);gb=f(b);while(b-a)>eps)c=(a+b)/2;gc=f(c);if(gc=0)break;else if(gc*gb<0)a=c;ga=gc;elseb=c;gb=gc;printx();printf("ttt方程的根为:X=%lfn",b);printx();return 0;9.霍纳法多项式求值算法代码 #include<stdio.h>#include<stdlib.h>void printx()int i=0

34、;for(i=0;i<20;i+)printf(" ");for(i=0;i<40;i+)printf("*");for(i=0;i<20;i+)printf(" ");main()int i=0,n=0;double *A=NULL,*B=NULL,c=0.0;printf("nntttt霍纳法求解多项式n");printx();printf("ttt请输入多项式的项数:");scanf("%d",&n);A=(double*)malloc(n*sizeof(double);B=(double*)malloc(n*sizeof(double);printf("nttt请输入各项的系数，按照由高次到低次排序:nttt");for(i=n-1;i>=0;i-)scanf("%lf",&Ai);printf("ttt请输入X的值:");scanf("%lf",&c);Bn-1=A