計算方法程序(C++)

1、1.經(jīng)典四階龍格庫塔法解一階微分方程組代碼#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;/微分方程的函數(shù)文件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為矩陣的行數(shù)。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);/下面是回代函數(shù)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請輸入待求解方程組的增廣矩陣的行數(shù):");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 /用來設置方程組的行數(shù)#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維數(shù)相同的列向量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(xiàn)(P)的誤差限，最大迭代次數(shù)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迭代次數(shù)為:%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("求解函數(shù)y=x*x在(0,1)上的龍貝格矩陣n");printf("龍貝格矩陣最大行數(shù)為%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("求解經(jīng)過點（0，0.0），（1，0.5），（2，2.0）和（3，1.5）n而且一階導數(shù)邊界條件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在區(qū)間(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請先輸入待擬合的點的個數(shù):");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請輸入需要擬合的次數(shù):");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次多項式系數(shù)為，冪次由高到低：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("迭代次數(shù)超過允許的最大

31、迭代次數(shù)！");printf("nnnttt不動點的近似值為:%lf",p);printf("nttt程序迭代次數(shù)為%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求根區(qū)間為（-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請輸入多項式的項數(shù):");scanf("%d",&n);A=(double*)malloc(n*sizeof(double);B=(double*)malloc(n*sizeof(double);printf("nttt請輸入各項的系數(shù)，按照由高次到低次排序:nttt");for(i=n-1;i>=0;i-)scanf("%lf",&Ai);printf("ttt請輸入X的值:");scanf("%lf",&c);Bn-1=A