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, long i);
  7. void FUNC(double x, double *y, double *f);
  8. void rk4m(double *y, double *f, int N, double a, double b, int n, void (*F)());
  9.  
  10. int main(void) {
  11.     int n = 100, N = 2;
  12.     double *y, *f, a = 0.0, b = 10.0;
  13.     y = dvector(1, N);
  14.     f = dvector(1, N);
  15.     y[1] = 0.1;
  16.     y[2] = 0.0;
  17.     rk4m(y, f, N, a, b, n, FUNC);
  18.     free_dvector(y, 1);
  19.     free_dvector(f, 1);
  20.     return 0;
  21. }
  22.  
  23. void rk4m(double *y, double *f, int N, double a, double b, int n, void (*F)()) {
  24.     double *k1, *k2, *k3, *k4, *tmp, x, h;
  25.     int i, j;
  26.     k1 = dvector(1, N);
  27.     k2 = dvector(1, N);
  28.     k3 = dvector(1, N);
  29.     k4 = dvector(1, N);
  30.     tmp = dvector(1, N);
  31.     h = (b - a) / n;
  32.     x = a;
  33.     printf("%.3lf\t %.3lf\t %.10lf\n", x, y[1], y[2]);
  34.     for (i = 0; i < n; i++) {
  35.         (*F)(x, y, f);
  36.         for (j = 1; j <= N; j++) k1[j] = h * f[j];
  37.         for (j = 1; j <= N; j++) tmp[j] = y[j] + k1[j] / 2.0;
  38.         (*F)(x + h / 2.0, tmp, f);
  39.         for (j = 1; j <= N; j++) k2[j] = h * f[j];
  40.         for (j = 1; j <= N; j++) tmp[j] = y[j] + k2[j] / 2.0;
  41.         (*F)(x + h / 2.0, tmp, f);
  42.         for (j = 1; j <= N; j++) k3[j] = h * f[j];
  43.         for (j = 1; j <= N; j++) tmp[j] = y[j] + k3[j];
  44.         (*F)(x + h, tmp, f);
  45.         for (j = 1; j <= N; j++) k4[j] = h * f[j];
  46.         for (j = 1; j <= N; j++) y[j] = y[j] + (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6.0;
  47.         x += h;
  48.         printf("%.3lf\t %.10lf\t %.10lf\n", x, y[1], y[2]);
  49.     }
  50.     free_dvector(k1, 1);
  51.     free_dvector(k2, 1);
  52.     free_dvector(k3, 1);
  53.     free_dvector(k4, 1);
  54.     free_dvector(tmp, 1);
  55. }
  56.  
  57. void FUNC(double x, double *y, double *f) {
  58.     f[1] = y[2];
  59.     f[2] = cos(x) - 4 * y[2] - 3 * y[1];
  60. }
  61.  
  62. double *dvector(long i, long j) {
  63.     double *a;
  64.     if ((a = malloc((j - i + 1) * sizeof(double))) == NULL) {
  65.         printf("メモリが確保できません\n");
  66.         exit(1);
  67.     }
  68.     return a - i;
  69. }
  70.  
  71. void free_dvector(double *a, long i) {
  72.     free((void *)(a + i));
  73. }
  74.  
Success #stdin #stdout 0.01s 5284KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
0.000	 0.100	 0.0000000000
0.100	 0.1030675008	 0.0572487294
0.200	 0.1107521521	 0.0933652025
0.300	 0.1212169032	 0.1136137864
0.400	 0.1330813637	 0.1219313032
0.500	 0.1453073605	 0.1212755212
0.600	 0.1571145645	 0.1138826631
0.700	 0.1679183270	 0.1014574965
0.800	 0.1772839019	 0.0853134699
0.900	 0.1848927370	 0.0664758373
1.000	 0.1905176381	 0.0457573632
1.100	 0.1940044361	 0.0238137211
1.200	 0.1952584048	 0.0011838597
1.300	 0.1942341310	 -0.0216807519
1.400	 0.1909278773	 -0.0443950935
1.500	 0.1853717292	 -0.0666256232
1.600	 0.1776290050	 -0.0880797634
1.700	 0.1677905456	 -0.1084986511
1.800	 0.1559716033	 -0.1276522644
1.900	 0.1423091277	 -0.1453362648
2.000	 0.1269593034	 -0.1613700587
2.100	 0.1100952349	 -0.1755957117
2.200	 0.0919047111	 -0.1878774357
2.300	 0.0725879994	 -0.1981014445
2.400	 0.0523556442	 -0.2061760200
2.500	 0.0314262518	 -0.2120316745
2.600	 0.0100242591	 -0.2156213198
2.700	 -0.0116223131	 -0.2169203815
2.800	 -0.0332841130	 -0.2159268071
2.900	 -0.0547327042	 -0.2126609382
3.000	 -0.0757427930	 -0.2071652187
3.100	 -0.0960944160	 -0.1995037280
3.200	 -0.1155750729	 -0.1897615262
3.300	 -0.1339817842	 -0.1780438123
3.400	 -0.1511230548	 -0.1644748933
3.500	 -0.1668207260	 -0.1491969723
3.600	 -0.1809116971	 -0.1323687618
3.700	 -0.1932495003	 -0.1141639354
3.800	 -0.2037057133	 -0.0947694303
3.900	 -0.2121711950	 -0.0743836171
4.000	 -0.2185571326	 -0.0532143542
4.100	 -0.2227958891	 -0.0314769454
4.200	 -0.2248416422	 -0.0093920216
4.300	 -0.2246708092	 0.0128166336
4.400	 -0.2222822518	 0.0349242908
4.500	 -0.2176972597	 0.0567074949
4.600	 -0.2109593127	 0.0779462730
4.700	 -0.2021336236	 0.0984263107
4.800	 -0.1913064651	 0.1179410732
4.900	 -0.1785842895	 0.1362938501
5.000	 -0.1640926479	 0.1532997047
5.100	 -0.1479749196	 0.1687873060
5.200	 -0.1303908660	 0.1826006270
5.300	 -0.1115150215	 0.1946004909
5.400	 -0.0915349378	 0.2046659506
5.500	 -0.0706492994	 0.2126954861
5.600	 -0.0490659295	 0.2186080102
5.700	 -0.0269997043	 0.2223436695
5.800	 -0.0046703986	 0.2238644352
5.900	 0.0176995174	 0.2231544757
6.000	 0.0398871069	 0.2202203088
6.100	 0.0616712003	 0.2150907304
6.200	 0.0828346098	 0.2078165219
6.300	 0.1031663045	 0.1984699381
6.400	 0.1224635229	 0.1871439807
6.500	 0.1405338033	 0.1739514655
6.600	 0.1571969095	 0.1590238915
6.700	 0.1722866354	 0.1425101241
6.800	 0.1856524684	 0.1245749044
6.900	 0.1971610957	 0.1053972009
7.000	 0.2066977390	 0.0851684190
7.100	 0.2141673031	 0.0640904861
7.200	 0.2194953280	 0.0423738325
7.300	 0.2226287350	 0.0202352868
7.400	 0.2235363581	 -0.0021040922
7.500	 0.2222092574	 -0.0244212253
7.600	 0.2186608090	 -0.0464932434
7.700	 0.2129265732	 -0.0680997154
7.800	 0.2050639399	 -0.0890248520
7.900	 0.1951515559	 -0.1090596622
8.000	 0.1832885406	 -0.1280040428
8.100	 0.1695934959	 -0.1456687785
8.200	 0.1542033219	 -0.1618774327
8.300	 0.1372718500	 -0.1764681118
8.400	 0.1189683062	 -0.1892950829
8.500	 0.0994756206	 -0.2002302303
8.600	 0.0789886005	 -0.2091643366
8.700	 0.0577119843	 -0.2160081736
8.800	 0.0358583957	 -0.2206933953
8.900	 0.0136462205	 -0.2231732200
9.000	 -0.0087025760	 -0.2234228989
9.100	 -0.0309646660	 -0.2214399633
9.200	 -0.0529175907	 -0.2172442493
9.300	 -0.0743419824	 -0.2108777006
9.400	 -0.0950237565	 -0.2024039488
9.500	 -0.1147562500	 -0.1919076781
9.600	 -0.1333422867	 -0.1794937797
9.700	 -0.1505961468	 -0.1652863033
9.800	 -0.1663454226	 -0.1494272182
9.900	 -0.1804327408	 -0.1320749948
10.000	 -0.1927173350	 -0.1134030214
https://ideone.com/mHeBIN
language:
C (gcc 8.3)
created:
3 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!