fork download
copy 
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. double *dvector(long i, long j);
  6. void free_dvector(double *a);
  7. void FUNC(double x, double *y, double *f);
  8. void rk4m(double *y, double *f, int N, double a, double b, double step, void (*F)(double, double *, double *));
  9.  
  10. int main(void) {
  11.     int N = 2;
  12.     double *y, *f, a = 0.0, b = 1.0, step = 0.01; // 終了値を 1.0、刻み幅を 0.01 に設定
  13.  
  14.     y = dvector(0, N - 1);
  15.     f = dvector(0, N - 1);
  16.  
  17.     y[0] = 0.1; 
  18.     y[1] = 0.0;
  19.  
  20.     rk4m(y, f, N, a, b, step, FUNC);
  21.  
  22.     free_dvector(y);
  23.     free_dvector(f);
  24.  
  25.     return 0;
  26. }
  27.  
  28. void FUNC(double x, double *y, double *f) {
  29.     f[0] = y[1]; // dy1/dx = y2
  30.     f[1] = cos(x) - 4 * y[1] - 3 * y[0]; // dy2/dx = cos(x) - 4y2 - 3y1
  31. }
  32.  
  33. void rk4m(double *y, double *f, int N, double a, double b, double step, void (*F)(double, double *, double *)) {
  34.     double *k1, *k2, *k3, *k4, *tmp, x, h;
  35.     int j;
  36.  
  37.     k1 = dvector(0, N - 1);
  38.     k2 = dvector(0, N - 1);
  39.     k3 = dvector(0, N - 1);
  40.     k4 = dvector(0, N - 1);
  41.     tmp = dvector(0, N - 1);
  42.  
  43.     h = step;
  44.     x = a;
  45.  
  46.     while (x <= b) {
  47.         printf("%.2lf\t %.10lf\t %.10lf\n", x, y[0], y[1]);
  48.  
  49.         F(x, y, f);
  50.         for (j = 0; j < N; j++) k1[j] = f[j];
  51.  
  52.         for (j = 0; j < N; j++) tmp[j] = y[j] + h * k1[j] / 2.0;
  53.         F(x + h / 2.0, tmp, f);
  54.         for (j = 0; j < N; j++) k2[j] = f[j];
  55.  
  56.         for (j = 0; j < N; j++) tmp[j] = y[j] + h * k2[j] / 2.0;
  57.         F(x + h / 2.0, tmp, f);
  58.         for (j = 0; j < N; j++) k3[j] = f[j];
  59.  
  60.         for (j = 0; j < N; j++) tmp[j] = y[j] + h * k3[j];
  61.         F(x + h, tmp, f);
  62.         for (j = 0; j < N; j++) k4[j] = f[j];
  63.  
  64.         for (j = 0; j < N; j++)
  65.             y[j] = y[j] + h * (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0;
  66.  
  67.         x += h; // x を更新
  68.     }
  69.  
  70.     free_dvector(k1);
  71.     free_dvector(k2);
  72.     free_dvector(k3);
  73.     free_dvector(k4);
  74.     free_dvector(tmp);
  75. }
  76.  
  77. double *dvector(long i, long j) {
  78.     double *a;
  79.     if ((a = malloc((j - i + 1) * sizeof(double))) == NULL) {
  80.         printf("メモリが確保できません\n");
  81.         exit(1);
  82.     }
  83.     return a;
  84. }
  85.  
  86. void free_dvector(double *a) {
  87.     free(a);
  88. }
  89.  
  90.  
Success #stdin #stdout 0.01s 5288KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
0.00	 0.1000000000	 0.0000000000
0.01	 0.1000345367	 0.0068613400
0.02	 0.1001363201	 0.0134506418
0.03	 0.1003026687	 0.0197756543
0.04	 0.1005309773	 0.0258439005
0.05	 0.1008187152	 0.0316626843
0.06	 0.1011634232	 0.0372390969
0.07	 0.1015627122	 0.0425800231
0.08	 0.1020142611	 0.0476921477
0.09	 0.1025158142	 0.0525819614
0.10	 0.1030651803	 0.0572557664
0.11	 0.1036602300	 0.0617196824
0.12	 0.1042988942	 0.0659796516
0.13	 0.1049791625	 0.0700414445
0.14	 0.1056990812	 0.0739106647
0.15	 0.1064567521	 0.0775927537
0.16	 0.1072503302	 0.0810929966
0.17	 0.1080780230	 0.0844165257
0.18	 0.1089380883	 0.0875683258
0.19	 0.1098288333	 0.0905532385
0.20	 0.1107486125	 0.0933759661
0.21	 0.1116958272	 0.0960410764
0.22	 0.1126689235	 0.0985530061
0.23	 0.1136663912	 0.1009160652
0.24	 0.1146867626	 0.1031344407
0.25	 0.1157286113	 0.1052121999
0.26	 0.1167905511	 0.1071532947
0.27	 0.1178712345	 0.1089615642
0.28	 0.1189693521	 0.1106407389
0.29	 0.1200836311	 0.1121944436
0.30	 0.1212128345	 0.1136262003
0.31	 0.1223557600	 0.1149394318
0.32	 0.1235112392	 0.1161374644
0.33	 0.1246781361	 0.1172235309
0.34	 0.1258553471	 0.1182007731
0.35	 0.1270417991	 0.1190722450
0.36	 0.1282364493	 0.1198409153
0.37	 0.1294382843	 0.1205096696
0.38	 0.1306463190	 0.1210813135
0.39	 0.1318595960	 0.1215585745
0.40	 0.1330771848	 0.1219441046
0.41	 0.1342981810	 0.1222404824
0.42	 0.1355217057	 0.1224502157
0.43	 0.1367469046	 0.1225757432
0.44	 0.1379729478	 0.1226194366
0.45	 0.1391990283	 0.1225836030
0.46	 0.1404243622	 0.1224704864
0.47	 0.1416481877	 0.1222822699
0.48	 0.1428697644	 0.1220210775
0.49	 0.1440883729	 0.1216889757
0.50	 0.1453033143	 0.1212879754
0.51	 0.1465139093	 0.1208200335
0.52	 0.1477194982	 0.1202870546
0.53	 0.1489194398	 0.1196908924
0.54	 0.1501131115	 0.1190333516
0.55	 0.1512999081	 0.1183161890
0.56	 0.1524792422	 0.1175411150
0.57	 0.1536505430	 0.1167097953
0.58	 0.1548132561	 0.1158238518
0.59	 0.1559668432	 0.1148848644
0.60	 0.1571107817	 0.1138943718
0.61	 0.1582445640	 0.1128538732
0.62	 0.1593676974	 0.1117648290
0.63	 0.1604797035	 0.1106286624
0.64	 0.1615801182	 0.1094467603
0.65	 0.1626684908	 0.1082204744
0.66	 0.1637443841	 0.1069511224
0.67	 0.1648073740	 0.1056399889
0.68	 0.1658570488	 0.1042883264
0.69	 0.1668930095	 0.1028973563
0.70	 0.1679148689	 0.1014682697
0.71	 0.1689222517	 0.1000022288
0.72	 0.1699147941	 0.0985003670
0.73	 0.1708921434	 0.0969637905
0.74	 0.1718539578	 0.0953935786
0.75	 0.1727999063	 0.0937907848
0.76	 0.1737296683	 0.0921564375
0.77	 0.1746429333	 0.0904915406
0.78	 0.1755394006	 0.0887970745
0.79	 0.1764187794	 0.0870739967
0.80	 0.1772807883	 0.0853232424
0.81	 0.1781251551	 0.0835457252
0.82	 0.1789516166	 0.0817423379
0.83	 0.1797599185	 0.0799139528
0.84	 0.1805498152	 0.0780614227
0.85	 0.1813210693	 0.0761855813
0.86	 0.1820734518	 0.0742872436
0.87	 0.1828067418	 0.0723672068
0.88	 0.1835207262	 0.0704262505
0.89	 0.1842151996	 0.0684651375
0.90	 0.1848899643	 0.0664846141
0.91	 0.1855448297	 0.0644854107
0.92	 0.1861796126	 0.0624682422
0.93	 0.1867941370	 0.0604338085
0.94	 0.1873882335	 0.0583827950
0.95	 0.1879617399	 0.0563158728
0.96	 0.1885145002	 0.0542336995
0.97	 0.1890463652	 0.0521369191
0.98	 0.1895571920	 0.0500261629
0.99	 0.1900468439	 0.0479020496
https://ideone.com/gTAMBk
微分方程式
language:
C (gcc 8.3)
created:
10 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!