01: #include <stdio.h>
02: #include <math.h>
03: #define NSTEPS 10000                // 計算ステップ数
04: #define NDIM 30000               // 用意する配列の大きさ
05: #define t_MAX 1.0                // 計算する最大 x
06: #define R 100
07: #define C 1e-6
08: #define V 1
09: #define omega 2*M_PI*50
10: 
11: double f(double t, double I);    // プロトタイプ宣言
12: double sol(double t);
13: 
14: //============================================================
15: // main 関数
16: //============================================================
17: int main(void){
18:   double t[NDIM], Sol_I[NDIM];  // 計算結果を入れる
19:   double Euler_I[NDIM], RK2_I[NDIM], RK4_I[NDIM];
20:   double dt, k1, k2, k3, k4;
21:   int i;
22:   FILE *result, *error;         // 計算結果を保存するファイル
23:   
24: 
25:   dt = t_MAX/NSTEPS;            // 計算のきざみ幅dxの計算
26:   t[0] = 0;
27:   Sol_I[0] = 0;
28:   Euler_I[0] = 0;
29:   RK2_I[0] = 0;
30:   RK4_I[0] = 0;
31: 
32:   result = fopen("result.txt","w"); // 結果を書くファイルをオープン
33:   error = fopen("error.txt","w"); // 結果を書くファイルをオープン
34: 
35:   fprintf(result,"%e\t%e\t%e\t%e\t%e\n", t[0],Sol_I[0],
36:       Euler_I[0],RK2_I[0],RK4_I[0]);
37:     
38:   //------- 微分方程式の計算 --------
39: 
40:   for(i=0; i<NSTEPS; i++){
41:     t[i+1] = t[0]+(i+1)*dt;       //x[i+1]=x[i]+dxよりも精度が良い
42:     Sol_I[i+1]= sol(t[i+1]);
43: 
44:     Euler_I[i+1] = Euler_I[i]+f(t[i], Euler_I[i])*dt;
45: 
46:     k1=f(t[i], RK2_I[i])*dt;
47:     k2=f(t[i]+dt, RK2_I[i]+k1)*dt;
48:     RK2_I[i+1] = RK2_I[i]+0.5*(k1+k2);
49: 
50:     k1=f(t[i], RK4_I[i])*dt;
51:     k2=f(t[i]+0.5*dt, RK4_I[i]+0.5*k1)*dt;
52:     k3=f(t[i]+0.5*dt, RK4_I[i]+0.5*k2)*dt;
53:     k4=f(t[i]+dt, RK4_I[i]+k3)*dt;
54:     RK4_I[i+1] = RK4_I[i]+(k1+2*k2+2*k3+k4)/6.0;
55: 
56:     fprintf(result,"%e\t%e\t%e\t%e\t%e\n", t[i+1],Sol_I[i+1],
57:         Euler_I[i+1],RK2_I[i+1],RK4_I[i+1]);
58: 
59:     fprintf(error,"%e\t%e\t%e\t%e\n", t[i+1],
60:         fabs(Euler_I[i+1]-Sol_I[i+1]),
61:         fabs(RK2_I[i+1]-Sol_I[i+1]),
62:         fabs(RK4_I[i+1]-Sol_I[i+1])
63:         );
64:   }
65:       
66:   fclose(result);
67:   fclose(error);
68:     
69:   return 0;
70: }
71: 
72: //============================================================
73: // dI/dt = f(I,t) の　f(I,t)を計算する関数
74: //============================================================
75: double f(double t, double I){
76:   return -I/(C*R)+omega*V/R*cos(omega*t);
77: }
78: 
79: //============================================================
80: // I(t)の解析解
81: //============================================================
82: double sol(double t){
83:   return V*omega*C*(cos(omega*t)+R*omega*C*sin(omega*t)-exp(-t/(C*R)))/
84:     (1+C*C*R*R*omega*omega);
85: }