欢迎来到天天文库
浏览记录
ID:38807346
大小:17.22 KB
页数:5页
时间:2019-06-19
《计算方法 上机操作(C语言)》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、#include"stdlib.h"#include"math.h"#include"stdio.h"#defineN4intagaus(doublea[N][N],doubleb[N],intn);//高斯全主元法求deltaX,具体见程序中注释voidGetFun(doublex[N],doubleb[N]);//F(X),结果放在b里voidGetJacobi(doublex[N],doubleJ[N][N]);//计算F'(X),结果放在J里doubleGetDeltaX(doubleJ[N][N],doubleb[N],doubledeltaX[N]);//求delta
2、X主程序,调用agausvoidUpdateX(doublex[N],doubledeltaX[N]);//xk+1=xk+deltaXvoidSolveNonlinear(doublex[N],doubleeps,intk);//非线性方程组主程序,迭代终止条件是最大迭代次数或者绝对值范数(街区距离)main(){inti,k;doubleeps;doublex[4]={1,1.0,1.0,1.0};doubleb[N];eps=0.0000001;k=100;SolveNonlinear(x,eps,k);GetFun(x,b);printf("");printf("
3、");for(i=0;i4、X);if((nIter>k)5、6、(NormDeltaX7、}voidGetJacobi(x,J)//计算F'(X),结果放在J里doublex[N],J[N][N];{J[0][0]=2*x[0];J[0][1]=2*x[1];J[0][2]=2*x[2];J[0][3]=2*x[3];J[1][0]=4*x[0];J[1][1]=2*x[1];J[1][2]=-4;J[1][3]=4*x[3];J[2][0]=6*x[0];J[2][1]=-4;J[2][2]=2*x[2];J[2][3]=6*x[3];J[3][0]=4;J[3][1]=2*x[1];J[3][2]=4*x[2];J[3][3]=10*x[3];return;}in8、tagaus(a,b,n)//高斯全主元法解方程,结果放在b里intn;doublea[N][N],b[N];{int*js,l,k,j,i,is;doubled,t;js=malloc(n*sizeof(int));//记载列交换标号l=1;for(k=0;k<=n-2;k++){d=0.0;for(i=k;i<=n-1;i++)//求全主元{for(j=k;j<=n-1;j++){t=fabs(a[i][j]);if(t>d){d=t;js[k]=j;//记住列is=i;//记住行}}}if(d+1.0==1.0)l=0;else{/*if(js[k]!=k)//列交换{fo9、r(i=0;i<=n-1;i++){t=a[i][k];a[i][k]=a[i][js[k]];a[i][js[k]]=t;}}*/if(is!=k){for(j=k;j<=n-1;j++)//行交换{t=a[k][j];a[k][j]=a[is][j];a[is][j]=t;}t=b[k];b[k]=b[is];b[is]=t;//b交换}}if(l==0)//最大主元为零不能计算出正确结果{free(js);printf("fail");return(0);}//以下为消
4、X);if((nIter>k)
5、
6、(NormDeltaX7、}voidGetJacobi(x,J)//计算F'(X),结果放在J里doublex[N],J[N][N];{J[0][0]=2*x[0];J[0][1]=2*x[1];J[0][2]=2*x[2];J[0][3]=2*x[3];J[1][0]=4*x[0];J[1][1]=2*x[1];J[1][2]=-4;J[1][3]=4*x[3];J[2][0]=6*x[0];J[2][1]=-4;J[2][2]=2*x[2];J[2][3]=6*x[3];J[3][0]=4;J[3][1]=2*x[1];J[3][2]=4*x[2];J[3][3]=10*x[3];return;}in8、tagaus(a,b,n)//高斯全主元法解方程,结果放在b里intn;doublea[N][N],b[N];{int*js,l,k,j,i,is;doubled,t;js=malloc(n*sizeof(int));//记载列交换标号l=1;for(k=0;k<=n-2;k++){d=0.0;for(i=k;i<=n-1;i++)//求全主元{for(j=k;j<=n-1;j++){t=fabs(a[i][j]);if(t>d){d=t;js[k]=j;//记住列is=i;//记住行}}}if(d+1.0==1.0)l=0;else{/*if(js[k]!=k)//列交换{fo9、r(i=0;i<=n-1;i++){t=a[i][k];a[i][k]=a[i][js[k]];a[i][js[k]]=t;}}*/if(is!=k){for(j=k;j<=n-1;j++)//行交换{t=a[k][j];a[k][j]=a[is][j];a[is][j]=t;}t=b[k];b[k]=b[is];b[is]=t;//b交换}}if(l==0)//最大主元为零不能计算出正确结果{free(js);printf("fail");return(0);}//以下为消
7、}voidGetJacobi(x,J)//计算F'(X),结果放在J里doublex[N],J[N][N];{J[0][0]=2*x[0];J[0][1]=2*x[1];J[0][2]=2*x[2];J[0][3]=2*x[3];J[1][0]=4*x[0];J[1][1]=2*x[1];J[1][2]=-4;J[1][3]=4*x[3];J[2][0]=6*x[0];J[2][1]=-4;J[2][2]=2*x[2];J[2][3]=6*x[3];J[3][0]=4;J[3][1]=2*x[1];J[3][2]=4*x[2];J[3][3]=10*x[3];return;}in
8、tagaus(a,b,n)//高斯全主元法解方程,结果放在b里intn;doublea[N][N],b[N];{int*js,l,k,j,i,is;doubled,t;js=malloc(n*sizeof(int));//记载列交换标号l=1;for(k=0;k<=n-2;k++){d=0.0;for(i=k;i<=n-1;i++)//求全主元{for(j=k;j<=n-1;j++){t=fabs(a[i][j]);if(t>d){d=t;js[k]=j;//记住列is=i;//记住行}}}if(d+1.0==1.0)l=0;else{/*if(js[k]!=k)//列交换{fo
9、r(i=0;i<=n-1;i++){t=a[i][k];a[i][k]=a[i][js[k]];a[i][js[k]]=t;}}*/if(is!=k){for(j=k;j<=n-1;j++)//行交换{t=a[k][j];a[k][j]=a[is][j];a[is][j]=t;}t=b[k];b[k]=b[is];b[is]=t;//b交换}}if(l==0)//最大主元为零不能计算出正确结果{free(js);printf("fail");return(0);}//以下为消
此文档下载收益归作者所有