《计算方法程序(C++)(精)》由会员分享，可在线阅读，更多相关《计算方法程序(C++)(精)（23页珍藏版）》请在金锄头文库上搜索。
1、1.经典四阶龙格库塔法解一阶微分方程组代码#include#include#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*)malloc(M+1)*sizeof(double);/T=0;Y=(double*
2、)malloc(M+1)*sizeof(double);/T=0;T0=a;for(i=1;iM+1;i+)Ti=Ti-1+h;Y0=ya;for(i=0;iM;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;iN+1;i+)for(j=0;j2;j+)printf(%lf ,Rji);printf(n);/return R;/微分方程的函数文件do
3、uble 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#include#include/在列向量中寻找绝对值最大的项，并返回该项的标号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)+
4、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
5、(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;printf(n增广矩阵高斯列主元消去后的矩阵为:n);for(j=0;j=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(
6、)int N=0,i=0;double *A=NULL,*X=NULL;printf(n请输入待求解方程组的增广矩阵的行数:);scanf(%d,&N);A=(double*)calloc(N*(N+1),sizeof(double);printf(请输入待求解方程组的增广矩阵(%d行%d列):n,N,N+1);for(i=0;iN*(N+1);i+)scanf(%lf,&Ai);system(cls);printf(方程的增广矩阵为:n);for(i=0;iN*(N+1);i+)if(i%(N+1)=0)printf(n);printf(%lft,Ai);uptrbk(A,N);X=back
7、sub(A,N);printf(nn方程组的解为:n);for(i=0;iN;i+)printf(X(%d)= %lfn,i+1,Xi);free(A);free(X);exit(0);3.牛顿法解非线性方程组算法程序代码#include#include#include#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(d
8、ouble 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=
9、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 *M
10、atrixMultiply(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=
11、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-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
