/* * testx-5.c --- 菊地文雄・山本昌宏『微分方程式と計算機演習』山海堂 * pp.125--127 のプログラムを C 言語に翻訳したものを * -(p u')'+r u=λu, u(0)=u(1)=0 * にした. */ #include #include #define NDIM (1000) #include #include int number_p = 0; int number_r = 0; double sqr(double x) { return x * x; } /* [0,1] 上の一様乱数 */ double drandom() { return random() / (double) INT_MAX; } /* プロトタイプ宣言 (定義はずっと後に出てくる) */ void TRID(double A[], double B[], double C[], double F[], int N, int ID); void INPR(double *X, double *Y, int N, double *S); void shinv(double a[], double b[], double c[], double x[], double y[], double eps, double h, double SHIFT, int N, int M, double *Eval, int *NITER, int *NPIVOT); void show_vector(int n, char *name, double x[]) { int i; for (i = 1; i <= n; i++) printf("%s[%d]=%f\n", name, i, x[i]); } double Pa, Pb, Ra, Rb, Rc, Oa, Ob, Oc, Od, Qa, Qb, Qc, Qd, Qe; double Ha, Hb, Ia, Ib, Ic, Ja, Jb, Jc, Jd; double P(double x){ if (number_p == 0) return 1.0; else if (number_p == 1) return Pa * x + Pb; else if (number_p == 2) return Ra * x * x + Rb * x + Rc; else if (number_p == 3) return Oa * x * x * x + Ob * x * x + Oc * x + Od; else if (number_p == 4) return Qa * x * x * x * x + Qb * x * x * x + Qc * x * x + Qd * x + Qe; } double R(double x){ if (number_r == 0) return 0.0; else if (number_r == 1) return Ha*x+Hb; else if (number_r == 2) return Ia * x * x + Ib * x + Ic; else if (number_r == 3) return Ja * x * x * x - Jb * x * x + Jc * x + Jd; } int main(void) { int N, M, i, k, di; double h, h2, h3, x, SHIFT, EPS, EVAL, MINEIGEN; int NITER, NPIVOT; double *A, *B, *C, *U, *V; char fname[128]; FILE *fp; double pi = 4 * atan(1.0); printf("INPUT : N=\n"); scanf("%d",&N); if(N<=1) exit(0); A = malloc((N+1) * sizeof(double)); B = malloc((N+1) * sizeof(double)); C = malloc((N+1) * sizeof(double)); U = malloc((N+1) * sizeof(double)); V = malloc((N+1) * sizeof(double)); EPS = 1e-15; M = 100; printf("0: P(x)=1.0,\n"); printf("1: P(x)=ax+b,\n"); printf("2: P(x)=ax^2+bx+c,\n"); printf("3: P(x)=ax^3+bx^2+cx+d, \n"); printf("4: P(x)=ax^4+bx^3+cx^2+dx+e \n"); printf("P(x)="); scanf("%d", &number_p); if (number_p == 1) { printf("a: "); scanf("%lf", &Pa); printf("b: "); scanf("%lf", &Pb); } else if (number_p == 2){ printf("a: "); scanf("%lf", &Ra); printf("b: "); scanf("%lf", &Rb); printf("c: "); scanf("%lf", &Rc); } else if (number_p == 3){ printf("a: "); scanf("%lf", &Oa); printf("b: "); scanf("%lf", &Ob); printf("c: "); scanf("%lf", &Oc); printf("d: "); scanf("%lf", &Od); } else if (number_p == 4){ printf("a: "); scanf("%lf", &Qa); printf("b: "); scanf("%lf", &Qb); printf("c: "); scanf("%lf", &Qc); printf("d: "); scanf("%lf", &Qd); printf("e: "); scanf("%lf", &Qe); } printf("0: R(x)=0,\n"); printf("1: R(x)=ax+b,\n"); printf("2: R(x)=ax^2+bx+c,\n"); printf("3: R(x)=ax^3+bx^2+cx+d \n"); printf("R(x)="); scanf("%d", &number_r); if (number_r == 1) { printf("a: "); scanf("%lf", &Ha); printf("b: "); scanf("%lf", &Hb); } else if (number_r == 2){ printf("a: "); scanf("%lf", &Ia); printf("b: "); scanf("%lf", &Ib); printf("c: "); scanf("%lf", &Ic); } else if (number_r == 3){ printf("a: "); scanf("%lf", &Ja); printf("b: "); scanf("%lf", &Jb); printf("c: "); scanf("%lf", &Jc); printf("d: "); scanf("%lf", &Jd); } printf("number p=%d\n", number_p); printf("number_r=%d\n", number_r); printf("n= "); scanf("%d", &di); h = 1.0/N; h2 = h*h; h3 = h/2; /* 計算結果をファイルに書き出す */ printf("file name: "); scanf("%s", fname); fp = fopen(fname, "w"); fprintf(fp, "# N, EPS, NITER\n"); fprintf(fp, "# %d %g %g %d\n", N, EPS, M); fprintf(fp, "# number_p=%d, number_r=%d\n",number_p,number_r); if (number_p == 9) { fprintf(fp, "# a=%f, b=%f\n", Pa, Pb); } for (k = 0; k < di; k++) { if (k == 0) { SHIFT = 0; } else { SHIFT = sqr((k + 1.0) / k) * EVAL; } for(i = 1; i < N; i++){ x = i * h; A[i] = -P(x-h3)/h2; B[i] = (P(x+h3)+P(x-h3)+h2*R(x))/h2-SHIFT; C[i] = -P(x+h3)/h2; } shinv(A,B,C,U,V,EPS,h,SHIFT,N-1,M,&EVAL,&NITER,&NPIVOT); printf("反復回数 = %3d EPS = %12.4e SHIFT = %12.4e NPIVOT = %3d\n", NITER, EPS, SHIFT, NPIVOT); if (k == 0) { MINEIGEN = EVAL; printf("固有値 = %20.15f\n",EVAL); fprintf(fp, "# 固有値 = %20.15f\n",EVAL); } else { printf("固有値 = %20.15f (%f)\n",EVAL, sqrt(EVAL / MINEIGEN)); fprintf(fp, "# 固有値 = %20.15f (%f)\n",EVAL, sqrt(EVAL / MINEIGEN)); } #ifdef OLD printf("** 節点での固有関数の値 **\n"); for(i=0; i<5; i++) printf(" I U(I) "); printf("\n"); for(i=1; i eps && *NITER < M); *NPIVOT=0; for (i=1; i<=N; i++) { if (b[i]<0.0) *NPIVOT=*NPIVOT+1; } } void TRID(double A[], double B[], double C[], double F[], int N, int ID) { int i; double AA; for(i=1; i< N; i++){ AA=A[i+1]/B[i]; if(ID==0) B[i+1]=B[i+1]-AA*C[i]; F[i+1]=F[i+1]-AA*F[i]; } F[N]=F[N]/B[N]; for(i=N-1; i>=1; i--) F[i]=(F[i]-C[i]*F[i+1])/B[i]; } void INPR(double *X, double *Y, int N, double *S) { int i; *S=0.0; for (i=1; i<=N; i++){ *S=*S+X[i]*Y[i]; } }