Metoda Rungego – Kutty

Mam do napisania program rozwiązujący, (a raczej "szukajacy" punktów) wykresu obliczanych za pomocą metody Rungego – Kutty. Ustaliłem sb krok na 0.2 ... i wiem tylko tyle. Wzory mam podane, ale z matmy jestem cienki, bo całek jeszcze nie mieliśmy (i różniczek, itp :/) i stoję w miejscu .. Mam to zrobić dla funkcji x(0) = 1, dx(t)/dt = exp(t). (e^x się znaczy).

Pomóżcie, proszę! Rozumiem tylko tyle, że moje punkty na osi x będę wyznaczał, poprzed dodanie kroku h,czyli:

x1 = 0, x2 = x1+h, x3 = x2+h, itd ... Co jednak z y?

Znalazłem taki kod:

/*
 
    Assume y(t) = exp(-2*t) then y(0) = 1.0 and dy/dt = -2*y(t) so we can use following
    FIRST-ORDER ODE for testing:
 
           dy/dt = exp(-2*t) - 3*y(t)
            y(0) = 1.0
 
*/
 
 
#include <stdio.h>
#include <math.h>
 
#define N     1            /* number of first order equations */
#define dist  0.005        /* stepsize in t*/
#define MAX   1.0          /* max for t */
 
FILE *output;              /* internal filename */
 
void runge4(double x, double y[], double step);   /* Runge-Kutta function */
double f(double x, double y[], int i);            /* function for derivatives */
 
main()
{
  double t, y[N], exactf,perc_diff;
  int j;
 
  output=fopen("output1.dat005", "w");            /* external filename */
 
  t   = 0.0;
  y[0]= 1.0;                                      /* initial position */
  exactf = exp (-2.0*t);
 
  fprintf(output, "%f\t%f\t%f\n", t, y[0],exactf);
 
  for (j=1; j*dist<=MAX ;j++)                     /* time loop */
     {
       t=j*dist;
       runge4(t, y, dist);
       exactf = exp (-2.0*t);
       perc_diff = (exactf-y[0])*100.0/exactf;
       fprintf(output, "%f\t%f\t%lf\t%lf\n", t, y[0],exactf,perc_diff);
     }
 
  fclose(output);
}
void runge4(double x, double y[], double step)
{
  double h=step/2.0,            /* the midpoint */
  t1[N], t2[N], t3[N],          /* temporary storage arrays */
  k1[N], k2[N], k3[N],k4[N];    /* for Runge-Kutta */
  int i;
 
  for (i=0;i<N;i++) t1[i]=y[i]+0.5*(k1[i]=step*f(x, y, i));
  for (i=0;i<N;i++) t2[i]=y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
  for (i=0;i<N;i++) t3[i]=y[i]+ (k3[i]=step*f(x+h, t2, i));
  for (i=0;i<N;i++) k4[i]= step*f(x+step, t3, i);
 
  for (i=0;i<N;i++) y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
 
double f(double x, double y[], int i)
{
  if (i==0) return(exp(-2*x) - 3*y[0]);     /* derivative of first equation */
  if (i!=0) exit(2);
}

no i niektóre części rozumiem .. te pętle, obliczają mi współrzędną y-kową dla funkcji exp(-2t), czy wystarczy, że zmienię w tym kodzie exp(-2t) na exp(t) i otrzymam to, co chciałbym otrzymać? :D

http://pl.wikipedia.org/wiki/Algorytm_Rungego-Kutty
Tam masz algorytm RK 4 rzędu
Masz parę (t,x) nie ma żadnego y.
Obliczasz po kolei współczynniki (na pcozatku za x(0) podstawiasz 1)
k1=f(tn,xn)=e^{tn (na początku będzie dla tn=0 e}0)
k2=f(tn+h/2,xn+1/2 hk1)= e^(tn+h/2)
k3...
k4...
Potem obliczasz dxn=1/6h(...)
i w końcu obliczasz
xn+1 = xn + dxn
To wszystko proste, dasz sam rade zaprogramować.

Ok, skoro nie ma żadnego y, to co będzie na OY na wykresie? O_o Nie bardzo rozumiem, jak to liczyć...

Coś takiego wymyśliłem:

1.pierwszy krok: (h = 0.01)

x(0) = 1
k1 = f(0, e^0)
k2 = f(0+1/2, e^0+1/20.01k1)
k3 = f(0+1/2, e^0+1/20.01k2)
k4 = f(0+h, e^0+0.01*k3)

y(0) = 1/6h(k1+2k2+2k3+k4)

i na wykresie rysuję punkt (x(0), y(0)), czyli (1, wynik_z_y(0)+x(0))

2.drugi krok:

x(1) = x(0)+h
k1 = f(x(0)+h, e^1)
k2 = f(x(0)+h+1/2, e^1+1/20.01k1)

... itd

czy to będzie ok? Ale jak obliczyć jeszcze [jak zakodzić, żeby poprawnie obliczało:]
coś takiego: f(x,y) ?

Weź mi powiedz, tak z ciekawości, gdzie jest tak, że najpierw się pisze programy dotyczące rachunku różniczkowego, a potem dopiero jest nauka tego od podstaw?
Co do twojego pytania, to nie jest kwestia różniczek, tylko nazewnictwa w funkcjach. Tłumacząc na tzw. chłopski rozum (bo teorię to w książkach znajdziesz):
y=f(x) to funkcja zmiennej x i tu masz te dwie osie, o których wspominasz. Literki t często używa się do oznaczenia czasu i wtedy funkcja zmiennej czasu, np. przebyta droga jest zapisywana
s=f(t) i masz dwie osie: pozioma t i pionowa s. Ty masz funkcję x=f(t), inaczej zapisane jako x(t), a więc t to oś pozioma, a x pionowa. To jest tylko kwestia nazewnictwa zmiennych.

ok, to jeszcze raz się upewnię: moje punkty będą wyglądać tak [napisze kilka przykładowych]: (0,e^{0), (1,e}1), (2, e^{2) itd... i te e}potega będą moimi y'. Tylko teraz, jak chcę obliczyć te k1,k2,k3,k4, to jak to będzie wyglądać? Czy to będzie tak:

0 krok: na wykresie rysuję punkt: (0,1)

1 krok:

k1 = e^0
k2 = e^(0+h/2)
k3 = e^(0+h/2)
k4 = e^(0+h)
y = 1/6h(k1+2k2+2k3+k4)

y(0+1) = y(1) = 1 + 1/6h(k1+2k2+2k3+k4)

i na wykresie rysuję punkt (0+h, 1 + 1/6h(k1+2k2+2k3+k4))

2 krok:

k1 = e^(0+h)
k2 = e^(0+h+h/2)
k3 = e^(0+h+h/2)
k4 = e^(0+h+h/2)
y = 1/6h(k1+2k2+2k3+k4)

y(1+1)=y(2) = [1 + 1/6h(k1+2k2+2k3+k4)] + 1/6h(k1+2k2+2k3+k4)]

i na wykresie rysuję punkt (0+h+h, 1 + 1/6h(k1+2k2+2k3+k4)] + 1/6h(k1+2k2+2k3+k4))

?

Wiem,że nikomu nie będzie chciało się spr wszystkiego, ale może mógłbyś zobaczyć, czy mniej więcej takie wyniki powinienem otrzymać?:

x są w porzadku, bo od 0.0 do 10.0:

 
y:  vector(1, 1.01005, 1.0202, 1.03045, 1.04081, 1.05127, 1.06184, 1.07251, 1.08329, 1.09417, 1.10517, 1.11628, 1.1275, 1.13883, 1.15027, 1.16183, 1.17351, 1.1853, 1.19722, 1.20925, 1.2214, 1.23368, 1.24608, 1.2586, 1.27125, 1.28403, 1.29693, 1.30996, 1.32313, 1.33643, 1.34986, 1.36343, 1.37713, 1.39097, 1.40495, 1.41907, 1.43333, 1.44773, 1.46228, 1.47698, 1.49182, 1.50682, 1.52196, 1.53726, 1.55271, 1.56831, 1.58407, 1.59999, 1.61607, 1.63232, 1.64872, 1.66529, 1.68203, 1.69893, 1.71601, 1.73325, 1.75067, 1.76827, 1.78604, 1.80399, 1.82212, 1.84043, 1.85893, 1.87761, 1.89648, 1.91554, 1.93479, 1.95424, 1.97388, 1.99372, 2.01375, 2.03399, 2.05443, 2.07508, 2.09594, 2.117, 2.13828, 2.15977, 2.18147, 2.2034, 2.22554, 2.24791, 2.2705, 2.29332, 2.31637, 2.33965, 2.36316, 2.38691, 2.4109, 2.43513, 2.4596, 2.48432, 2.50929, 2.53451, 2.55998, 2.58571, 2.6117, 2.63794, 2.66446, 2.69123, 2.71828, 2.7456, 2.77319, 2.80107, 2.82922, 2.85765, 2.88637, 2.91538, 2.94468, 2.97427, 3.00417, 3.03436, 3.06485, 3.09566, 3.12677, 3.15819, 3.18993, 3.22199, 3.25437, 3.28708, 3.32012, 3.35348, 3.38719, 3.42123, 3.45561, 3.49034, 3.52542, 3.56085, 3.59664, 3.63279, 3.6693, 3.70617, 3.74342, 3.78104, 3.81904, 3.85743, 3.89619, 3.93535, 3.9749, 4.01485, 4.0552, 4.09596, 4.13712, 4.1787, 4.2207, 4.26311, 4.30596, 4.34924, 4.39295, 4.4371, 4.48169, 4.52673, 4.57223, 4.61818, 4.66459, 4.71147, 4.75882, 4.80665, 4.85496, 4.90375, 4.95303, 5.00281, 5.05309, 5.10387, 5.15517, 5.20698, 5.25931, 5.31217, 5.36556, 5.41948, 5.47395, 5.52896, 5.58453, 5.64065, 5.69734, 5.7546, 5.81244, 5.87085, 5.92986, 5.98945, 6.04965, 6.11045, 6.17186, 6.23389, 6.29654, 6.35982, 6.42374, 6.4883, 6.5535, 6.61937, 6.68589, 6.75309, 6.82096, 6.88951, 6.95875, 7.02869, 7.09933, 7.17068, 7.24274, 7.31553, 7.38906, 7.46332, 7.53832, 7.61409, 7.69061, 7.7679, 7.84597, 7.92482, 8.00447, 8.08492, 8.16617, 8.24824, 8.33114, 8.41487, 8.49944, 8.58486, 8.67114, 8.75828, 8.84631, 8.93521, 9.02501, 9.11572, 9.20733, 9.29987, 9.39333, 9.48774, 9.58309, 9.6794, 9.77668, 9.87494, 9.97418, 10.0744, 10.1757, 10.2779, 10.3812, 10.4856, 10.591, 10.6974, 10.8049, 10.9135, 11.0232, 11.134, 11.2459, 11.3589, 11.473, 11.5883, 11.7048, 11.8224, 11.9413, 12.0613, 12.1825, 12.3049, 12.4286, 12.5535, 12.6797, 12.8071, 12.9358, 13.0658, 13.1971, 13.3298, 13.4637, 13.5991, 13.7357, 13.8738, 14.0132, 14.154, 14.2963, 14.44, 14.5851, 14.7317, 14.8797, 15.0293, 15.1803, 15.3329, 15.487, 15.6426, 15.7998, 15.9586, 16.119, 16.281, 16.4446, 16.6099, 16.7769, 16.9455, 17.1158, 17.2878, 17.4615, 17.637, 17.8143, 17.9933, 18.1741, 18.3568, 18.5413, 18.7276, 18.9158, 19.106, 19.298, 19.4919, 19.6878, 19.8857, 20.0855, 20.2874, 20.4913, 20.6972, 20.9052, 21.1153, 21.3276, 21.5419, 21.7584, 21.9771, 22.198, 22.421, 22.6464, 22.874, 23.1039, 23.3361, 23.5706, 23.8075, 24.0468, 24.2884, 24.5325, 24.7791, 25.0281, 25.2797, 25.5337, 25.7903, 26.0495, 26.3113, 26.5758, 26.8429, 27.1126, 27.3851, 27.6604, 27.9383, 28.2191, 28.5027, 28.7892, 29.0785, 29.3708, 29.666, 29.9641, 30.2652, 30.5694, 30.8766, 31.187, 31.5004, 31.817, 32.1367, 32.4597, 32.7859, 33.1155, 33.4483, 33.7844, 34.124, 34.4669, 34.8133, 35.1632, 35.5166, 35.8735, 36.2341, 36.5982, 36.9661, 37.3376, 37.7128, 38.0918, 38.4747, 38.8613, 39.2519, 39.6464, 40.0448, 40.4473, 40.8538, 41.2644, 41.6791, 42.098, 42.5211, 42.9484, 43.3801, 43.816, 44.2564, 44.7012, 45.1504, 45.6042, 46.0625, 46.5255, 46.9931, 47.4654, 47.9424, 48.4242, 48.9109, 49.4024, 49.899, 50.4004, 50.907, 51.4186, 51.9354, 52.4573, 52.9845, 53.517, 54.0549, 54.5982, 55.1469, 55.7011, 56.2609, 56.8263, 57.3975, 57.9743, 58.557, 59.1455, 59.7399, 60.3403, 60.9467, 61.5592, 62.1779, 62.8028, 63.434, 64.0715, 64.7155, 65.3659, 66.0228, 66.6863, 67.3565, 68.0335, 68.7172, 69.4079, 70.1054, 70.81, 71.5216, 72.2404, 72.9665, 73.6998, 74.4405, 75.1886, 75.9443, 76.7075, 77.4785, 78.2571, 79.0436, 79.838, 80.6404, 81.4509, 82.2695, 83.0963, 83.9314, 84.7749, 85.6269, 86.4875, 87.3567, 88.2347, 89.1214, 90.0171, 90.9218, 91.8356, 92.7586, 93.6908, 94.6324, 95.5835, 96.5441, 97.5144, 98.4944, 99.4843, 100.484, 101.494, 102.514, 103.544, 104.585, 105.636, 106.698, 107.77, 108.853, 109.947, 111.052, 112.168, 113.296, 114.434, 115.584, 116.746, 117.919, 119.104, 120.301, 121.51, 122.732, 123.965, 125.211, 126.469, 127.74, 129.024, 130.321, 131.631, 132.954, 134.29, 135.639, 137.003, 138.38, 139.77, 141.175, 142.594, 144.027, 145.474, 146.936, 148.413, 149.905, 151.411, 152.933, 154.47, 156.022, 157.591, 159.174, 160.774, 162.39, 164.022, 165.67, 167.335, 169.017, 170.716, 172.431, 174.164, 175.915, 177.683, 179.469, 181.272, 183.094, 184.934, 186.793, 188.67, 190.566, 192.481, 194.416, 196.37, 198.343, 200.337, 202.35, 204.384, 206.438, 208.513, 210.608, 212.725, 214.863, 217.022, 219.203, 221.406, 223.632, 225.879, 228.149, 230.442, 232.758, 235.097, 237.46, 239.847, 242.257, 244.692, 247.151, 249.635, 252.144, 254.678, 257.238, 259.823, 262.434, 265.072, 267.736, 270.426, 273.144, 275.889, 278.662, 281.463, 284.291, 287.149, 290.035, 292.949, 295.894, 298.867, 301.871, 304.905, 307.969, 311.064, 314.191, 317.348, 320.538, 323.759, 327.013, 330.3, 333.619, 336.972, 340.359, 343.779, 347.234, 350.724, 354.249, 357.809, 361.405, 365.037, 368.706, 372.412, 376.155, 379.935, 383.753, 387.61, 391.506, 395.44, 399.415, 403.429, 407.483, 411.579, 415.715, 419.893, 424.113, 428.375, 432.681, 437.029, 441.421, 445.858, 450.339, 454.865, 459.436, 464.054, 468.717, 473.428, 478.186, 482.992, 487.846, 492.749, 497.701, 502.703, 507.755, 512.859, 518.013, 523.219, 528.477, 533.789, 539.153, 544.572, 550.045, 555.573, 561.157, 566.796, 572.493, 578.246, 584.058, 589.928, 595.857, 601.845, 607.894, 614.003, 620.174, 626.407, 632.702, 639.061, 645.484, 651.971, 658.523, 665.142, 671.826, 678.578, 685.398, 692.287, 699.244, 706.272, 713.37, 720.539, 727.781, 735.095, 742.483, 749.945, 757.482, 765.095, 772.784, 780.551, 788.396, 796.319, 804.322, 812.406, 820.571, 828.818, 837.147, 845.561, 854.059, 862.642, 871.312, 880.069, 888.914, 897.847, 906.871, 915.985, 925.191, 934.489, 943.881, 953.367, 962.949, 972.626, 982.401, 992.275, 1002.25, 1012.32, 1022.49, 1032.77, 1043.15, 1053.63, 1064.22, 1074.92, 1085.72, 1096.63, 1107.65, 1118.79, 1130.03, 1141.39, 1152.86, 1164.45, 1176.15, 1187.97, 1199.91, 1211.97, 1224.15, 1236.45, 1248.88, 1261.43, 1274.11, 1286.91, 1299.84, 1312.91, 1326.1, 1339.43, 1352.89, 1366.49, 1380.22, 1394.09, 1408.1, 1422.26, 1436.55, 1450.99, 1465.57, 1480.3, 1495.18, 1510.2, 1525.38, 1540.71, 1556.2, 1571.84, 1587.63, 1603.59, 1619.71, 1635.98, 1652.43, 1669.03, 1685.81, 1702.75, 1719.86, 1737.15, 1754.61, 1772.24, 1790.05, 1808.04, 1826.21, 1844.57, 1863.11, 1881.83, 1900.74, 1919.85, 1939.14, 1958.63, 1978.31, 1998.2, 2018.28, 2038.56, 2059.05, 2079.74, 2100.65, 2121.76, 2143.08, 2164.62, 2186.37, 2208.35, 2230.54, 2252.96, 2275.6, 2298.47, 2321.57, 2344.9, 2368.47, 2392.27, 2416.32, 2440.6, 2465.13, 2489.91, 2514.93, 2540.2, 2565.73, 2591.52, 2617.57, 2643.87, 2670.44, 2697.28, 2724.39, 2751.77, 2779.43, 2807.36, 2835.57, 2864.07, 2892.86, 2921.93, 2951.3, 2980.96, 3010.92, 3041.18, 3071.74, 3102.61, 3133.79, 3165.29, 3197.1, 3229.23, 3261.69, 3294.47, 3327.58, 3361.02, 3394.8, 3428.92, 3463.38, 3498.19, 3533.34, 3568.85, 3604.72, 3640.95, 3677.54, 3714.5, 3751.83, 3789.54, 3827.63, 3866.09, 3904.95, 3944.19, 3983.83, 4023.87, 4064.31, 4105.16, 4146.42, 4188.09, 4230.18, 4272.69, 4315.64, 4359.01, 4402.82, 4447.07, 4491.76, 4536.9, 4582.5, 4628.55, 4675.07, 4722.06, 4769.52, 4817.45, 4865.87, 4914.77, 4964.16, 5014.05, 5064.45, 5115.34, 5166.75, 5218.68, 5271.13, 5324.11, 5377.61, 5431.66, 5486.25, 5541.39, 5597.08, 5653.33, 5710.15, 5767.53, 5825.5, 5884.05, 5943.18, 6002.91, 6063.24, 6124.18, 6185.73, 6247.9, 6310.69, 6374.11, 6438.17, 6502.88, 6568.23, 6634.24, 6700.92, 6768.26, 6836.29, 6904.99, 6974.39, 7044.48, 7115.28, 7186.79, 7259.02, 7331.97, 7405.66, 7480.09, 7555.27, 7631.2, 7707.89, 7785.36, 7863.6, 7942.63, 8022.46, 8103.08, 8184.52, 8266.78, 8349.86, 8433.78, 8518.54, 8604.15, 8690.62, 8777.97, 8866.19, 8955.29, 9045.29, 9136.2, 9228.02, 9320.77, 9414.44, 9509.06, 9604.62, 9701.15, 9798.65, 9897.13, 9996.6, 10097.1, 10198.5, 10301, 10404.6, 10509.1, 10614.8, 10721.4, 10829.2, 10938, 11047.9, 11159, 11271.1, 11384.4, 11498.8, 11614.4, 11731.1, 11849, 11968.1, 12088.4, 12209.9, 12332.6, 12456.5, 12581.7, 12708.2, 12835.9, 12964.9, 13095.2, 13226.8, 13359.7, 13494, 13629.6, 13766.6, 13904.9, 14044.7, 14185.8, 14328.4, 14472.4, 14617.9, 14764.8, 14913.2, 15063, 15214.4, 15367.3, 15521.8, 15677.8, 15835.3, 15994.5, 16155.2, 16317.6, 16481.6, 16647.2, 16814.6, 16983.5, 17154.2, 17326.6, 17500.8, 17676.7, 17854.3, 18033.7, 18215, 18398.1, 18583, 18769.7, 18958.4, 19148.9, 19341.3, 19535.7, 19732.1, 19930.4, 20130.7, 20333, 20537.3, 20743.7, 20952.2, 21162.8, 21375.5, 21590.3, 21807.3, 22026.5, 22247.8)