| //第12回 微分方程式の数値演算(1) 課題1 #include <stdio.h> #include <math.h> double func(double t, double x){ double a = 0.1; //[1/sec] double b = 0.02; //[1/m] double g = 9.807; //[m/sec2] return -b*x*x - a*x + g; } double rkm(double (*fnc)(double, double), double t, double x, double h){ //ルンゲクッタ法 double ck0, ck1, ck2, ck3; ck0=(*fnc)(t, x)*h; ck1=(*fnc)(t+h*0.5, x+ck0*0.5)*h; ck2=(*fnc)(t+h*0.5, x+ck1*0.5)*h; ck3=(*fnc)(t+h, x+ck2)*h; return (x+(ck0+(ck1+ck2)*2.0+ck3)/6.0); } double rkgm(double (*fnc)(double, double), double t, double x, double h, double *q){ //ルンゲクッタ・ジル法 static double p1=0.2928932188134524, p2=1.70710678118547; double ck0, ck1, ck2, ck3; double q1, q2, q3; double r1, r2, r3, r4; double tt, x1, x2, x3; ck0=(*fnc)(t,x)*h; r1=ck0*0.5-*q; x1=x+r1; q1=*q+r1*3.0-ck0*0.5; tt=t+h *0.5; ck1=(*fnc)(tt ,x1)*h; r2=(ck1-q1 )*p1; x2=x1+ r2; q2=q1+r2*3.0 -ck1*p1; ck2=(*fnc )(tt,x2)*h; r3=(ck2-q2 )*p2; x3=x2+r3; q3=q2+r3*3.0 -ck2*p2; tt=t+h; ck3=(*fnc)( tt,x3)*h; r4=(ck3-q3-q3 )/6.0; *q=q3 +r4*3.0-ck3*0.5; return (x3+r4); } double mlnm(double (*fnc) (double, double), double df[ ], double t, double x[], double h, double eps){ //ミルンの方法 double mlnm, ff , px; px=x[0]+((df[3]+df[1])*2.0-df[2])*h*4.0/3.0; while(1){ ff=(*fnc)(t, px); mlnm=x[2]+(ff+df[3]*4.0+df[2])*h/3.0; if(fabs(px-mlnm)< =eps){ return mlnm; } px=mlnm; } } double admsm(double (*fnc)(double, double), double df[ ], double t, double x[], double h, double eps){ //アダムス法 double admsm, ff, hd, px; hd= h/24.0; px=x[3]+(df[3]*55.0-df[2]*59.0+df [1]*37.0-df[0]*9.0)*hd; while(1){ ff=( *fnc)(t, px); admsm=x[3]+( ff*9.0+df[3]*19.0-df[2]*5.0+df[1])*hd; if(fabs(px-admsm) <=eps){ return admsm; } px = admsm ; } } void main(){ double eps = 0.00001, h = 0.1, xx = 10.0, q = 0; double ｘ[5] = {0}, df[5] = {0}; double (*f1)(double,double); double rk = 0, rkg = 0, mln = 0, adms = 0; int j = 0; double t = 0, t0 = 0; f1 = &func; printf("ルンゲクッタ法\nt,v\n" ); for(double i = 0; i <= xx; i += h){ printf("%f,%f\n", i, rk); rk = rkm(f1, i, rk, h); if(j <= 3){ x[j] = rk; df[j] = func( i, x[j]); } j++; } printf("\n"); printf("ルンゲクッタ・ギル法\nt,v\n"); for(double i = 0; i <= xx; i += h){ printf("%f,%f\n", i, rkg); rkg = rkgm(f1 , i, rkg, h, &q); } printf("\n"); printf("ミルン法\nt,v\n0,%f\n",mln); j = 0 ; for(double i = 0; i <= 3*h; i += h){ printf("%f,%f\n", i+h, x[j]); j++; } for(double i = 3*h; i < xx; i += h){ mln = mlnm(f1, df, i, x, h, eps); x[4] = mln; df[4] = func(i, mln); for(int k = 0; k <= 3; k++){ x[k] = x[k+1]; df[k] = df[k+1]; } x[0] = x[1]; x[1] = x[2]; x[2] = x[3]; x[3] = mln; df[0] = df[1]; df[1] = df[2]; df[2] = df[3]; df[3] = func(i, x[3]); printf("%f,%f\n", i+h, mln); } printf("\n"); printf("アダムス法\nt,v\n"); j = 0; adms = 0; x[0] = 0; for(int i=0; i<=3; i++){ t=t0+i*h; df[i]=func(t, x[i]); if(i<3) { x[i+1]=rkm(f1, t, x[i], h); } printf("%f,%f\n", t, x[i]); } for(double i = 3*h; i <= xx; i += h){ adms = admsm(f1, df, i, x, h, eps); x[4] = adms; df[4] = func(i, adms); for(int k = 0; k <= 3; k++){ x[k] = x[k+1]; df[k ] = df[k+1]; } printf("%f,%f\n", i , adms); } printf("\n"); } |
|
//第12回 微分方程式の数値演算(1) 課題2 |