buaa数值分析第二次大作业

#include<stdio.h> #include<math.h> #include<stdlib.h> #define N 10 #define E pow(10,-12) double A[10][10]; void putin(){ int i,j; for(i = 1 ; i <= N ; i++) { for(j = 1 ; j <= N ; j++){ if(i=j) A[i-1][j-1] = 1.52 * cos((double)(i + 1.2 * j)); else A[i-1][j-1] = sin((double)(0.5 * i + 0.2 * j)); } } for(i=0;i<N;i++){ for(j=0;j<N;j++) printf("%.12E ",A[i][j]); printf("\n"); } } void putoff(double A[N][N]){ int i,j; } int main(){ putin(); return 0 ; } /* double mifa(){ int i,k=0; double v,w; double bata[2000]={0}; double u[N],y[N],Y[N+4]={0}; for(i=0;i<N;i++) u[i]=1; // for(k=1;((bata[k]-bata[k-1])/bata[k]>E)||((bata[k]-bata[k-1])/bata[k]<(-E));k++) do{ v=0; k++; for(i=0;i<N;i++) v+=u[i]*u[i]; w=sqrt(v); for(i=0;i<N;i++) y[i]=u[i]/w; for(i=0;i<N;i++) Y[i+2]=y[i]; for(i=0;i<N;i++) u[i]=c*Y[i]+b*Y[i+1]+C[2][i]*Y[i+2]+b*Y[i+3]+c*Y[i+4]; for(i=0;i<N;i++) bata[k]+=y[i]*u[i]; // printf("bata(%d)= %.12E\n",k,bata[k]); }while((fabs(bata[k]-bata[k-1])/fabs(bata[k]))>E); // printf("%.12E\n",bata[k]); return bata[k]; } int max_2(int a, int b)//求两个数最大 { return(a>b ? a : b); } int min_2(int a, int b)//求两个数最小 { return(a<b ? a : b); } int max_3(int a, int b, int c)//求三个数最大 { int t; if (a>b) t = a; else t = b; return(t>c ? t : c); } void DLT() { int i, j, k, t; double p = 0; double q = 0; for (k = 0; k <N; k++) { for (j = k; (k!=0)&&(j < min_2(k + s+1, N)); j++) { p = 0; for (t = max_3(0, k - r, j - s); t < k ; t++) p += C[k - t + s][t ] * C[t - j + s][j ]; C[k-j+s][j] = C[k-j+s][j] - p; } for (i = k + 1; i < min_2(k + r + 1, N); i++) { q = 0; if(k==0) C[i-k+s][k]=C[i-k+s][k]/C[s][k]; else { for (t = max_3(0, i - r, k - s); t < k ; t++) q += C[i - t + s][t ] * C[t - k + s][k ]; C[i - k + s][k ] = (C[i - k + s][k ] - q) / C[s][k ]; } } } } double fmifa(){ int i , k = 0 , t , j; double v , w,sum,sum2; double bata[2000] = {0} , lambda[2000]={0}; double u[N] , y[N],Y[N]; for(i = 0 ; i < N ; i++) u[i] = 1; // for(k=1;((bata[k]-bata[k-1])/bata[k]>E)||((bata[k]-bata[k-1])/bata[k]<(-E));k++) DLT(); do{ v = 0; k++; for(i = 0;i < N ; i++) v+=u[i]*u[i]; w=sqrt(v); for(i=0;i<N;i++) y[i]=u[i]/w; for(i=0; i < N; i++) Y[i] = y[i]; //求 解 L U for(i = 1 ; i < N ; i++){ sum = 0; for( t = max_2 (0 , i-r ) ; t < i; t++) sum += C[i - t + s][ t]*y [t]; y[i]=y[i]-sum; } u[N-1]=y[N-1]/C[s][N-1]; for(j = N-2 ;j >= 0 ; j-- ){ sum2 = 0; for(t = j+1 ; t <= min_2(j+s,N-1) ; t++ ) sum2+=C[j-t+s][t]*u[t] ; u[j]= ( y[j] - sum2 ) / C[s][j]; } for(i = 0 ; i < N ; i++) bata[k] += Y[i]*u[i] ; lambda[k] = 1/bata[k]; // printf("lambda(%d)= %.12E\n",k,lambda[k]); }while((fabs(lambda[k]-lambda[k-1])/fabs(lambda[k]))>E); // printf("%.12E\n",bata[k]); return lambda[k]; } int main(){ int i,k; double cond,lambda_max,lambda1,lambda501,lambda_s,det = 1; double miu[39]={0},near_lambda[39]={0}; putin(); // putoff(C); lambda_max = mifa(); // printf("|lambda|_max = %.12E",lambda_max); if(lambda_max>=0){ lambda501=lambda_max; for(i=0;i<N;i++) C[2][i]=C[2][i]-lambda501; lambda1=mifa() + lambda501; } else{ lambda1=lambda_max; for(i=0;i<N;i++) C[2][i]=C[2][i]-lambda1; lambda501=mifa() + lambda1; } printf("lambda1=%.12E\n",lambda1); printf("lambda501=%.12E\n",lambda501); putin(); lambda_s=fmifa(); printf("lambda_s=%.12E\n",lambda_s); for( k = 1 ; k <= 39;k++){ miu[k] = lambda1+k*(lambda501 - lambda1) / 40 ; putin(); for( i =0;i < N ; i++) C[2][i]=C[2][i]-miu[k]; near_lambda[k] = fmifa() + miu[k]; printf("near_lambda_%d = %.12E\n",k,near_lambda[k]); } cond = max_2(fabs(lambda1/lambda_s), fabs(lambda501/lambda_s)); printf("cond(A)_2=%.12E\n",cond); putin(); DLT(); for (i = 0; i < N; i++) { det *= C[2][i]; } printf("det(A)=%.12E\n",det); return 0; } */