LCR過渡応答プログラム

---

プログラム

001: #include <stdio.h>
002: #include <math.h>
003: #define IMAX 100001
004: 
005: double f1(double I0, double I1);
006: double f2(double I0, double I1);
007: 
008: int mk_plot_data(char *a, double x[], double y[], int n);
009: int mk_graph(char *f);
010: 
011: double L=1.0;
012: double C=1.0;
013: double R=1.0;
014: 
015: /*================================================================*/
016: /*       main function                                            */
017: /*================================================================*/
018: int main(){
019:   double t[IMAX], I0[IMAX], I1[IMAX];
020:   double final_t, h;
021:   double k01, k02, k03, k04;
022:   double k11, k12, k13, k14;
023:   char temp;
024:   int ncal, i;
025: 
026: 
027:   /*--- set initial condition and cal range ---*/
028: 
029:   t[0] = 0.0;
030:   I0[0] = 0.0;
031:   I1[0] = 1.0;
032: 
033:   printf("\nfinal t = ");
034:   scanf("%lf%c", &final_t, &temp);
035: 
036:   do{
037:     printf("\nNumber of calculation steps( <= %d ) = ",IMAX-1);
038:     scanf("%d%c", &ncal, &temp);
039:   }while(ncal > IMAX-1);      
040: 
041: 
042:   /* --- size of calculation step --- */
043: 
044:   h=(final_t-t[0])/ncal;
045: 
046:   /* --- 4th Runge Kutta Calculation --- */
047: 
048:   for(i=0; i < ncal; i++){
049: 
050:     k01 = h*f1(I0[i],I1[i]);
051:     k11 = h*f2(I0[i],I1[i]);
052: 
053:     k02 = h*f1(I0[i]+k01/2, I1[i]+k11/2);
054:     k12 = h*f2(I0[i]+k01/2, I1[i]+k11/2);
055: 
056:     k03 = h*f1(I0[i]+k02/2, I1[i]+k12/2);
057:     k13 = h*f2(I0[i]+k02/2, I1[i]+k12/2);
058: 
059:     k04 = h*f1(I0[i]+k03, I1[i]+k13);
060:     k14 = h*f2(I0[i]+k03, I1[i]+k13);
061: 
062:     t[i+1] = t[i]+h;
063:     I0[i+1] = I0[i]+1.0/6.0*(k01+2.0*k02+2.0*k03+k04);
064:     I1[i+1] = I1[i]+1.0/6.0*(k11+2.0*k12+2.0*k13+k14);
065: 
066:   }
067: 
068: 
069:   /* --- make a graph --- */
070: 
071:   mk_plot_data("out.txt", t, I0, ncal);
072:   mk_graph("out.txt");
073: 
074:   return 0;
075: }
076: 
077: 
078: /*================================================================*/
079: /*       define function dI_0/dt                                  */
080: /*================================================================*/
081: double f1(double I0, double I1){
082:   double dI0dt;
083: 
084:   dI0dt = I1;
085: 
086:   return(dI0dt);
087: }
088: 
089: 
090: /*================================================================*/
091: /*       define function dI_1/dt                                  */
092: /*================================================================*/
093: double f2(double I0, double I1){
094:   double dI1dt;
095: 
096:   dI1dt = -1.0/L*(I0/C+R*I1);
097: 
098:   return(dI1dt);
099: }
100: 
101: 
102: /*==========================================================*/
103: /*   make a data file                                       */
104: /*==========================================================*/
105: int mk_plot_data(char *a, double x[], double y[], int n){
106:   int i;
107:   FILE *out;
108: 
109:   out = fopen(a, "w");
110: 
111:   for(i=0; i<=n; i++){
112:     fprintf(out, "%e\t%e\n", x[i], y[i]);
113:   }
114: 
115:   fclose(out);
116:   
117:   return 0;
118: }
119: 
120: /*==========================================================*/
121: /*   make a graph                                           */
122: /*==========================================================*/
123: int mk_graph(char *f){
124:   FILE *gp;
125:      
126:   gp = popen("gnuplot -persist","w");
127:       
128:   fprintf(gp, "reset\n");
129:   fprintf(gp, "set terminal postscript eps color\n");
130:   fprintf(gp, "set output \"graph.eps\"\n");
131:   fprintf(gp, "set nokey\n");
132:   fprintf(gp, "set grid\n");
133: 
134: 
135:   /* -------  plat graph ---------*/
136: 
137:   fprintf(gp, "plot \"%s\" using 1:2 with line\n",f);
138:   fprintf(gp, "set terminal x11\n");
139:   fprintf(gp, "replot\n");
140:                     
141:   pclose(gp);
142: 
143:   return 0;
144: }

---

last update: