RungeKutta.mw

Computing iterated derivatives of a function computed with a Runge - Kutta operator 

 

Clara Masse, John Masse and François Ollivier 

 

We define a Runge - Kutta operator for the differential equation y'=-y^2. 

> RK := proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*...
RK := proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*...
RK := proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*...
RK := proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*...
RK := proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*...
RK := proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*...
RK := proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*...
RK := proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*...
 

proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `...
proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `...
proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `...
proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `...
proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `...
proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `...
proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `...
proc (y0, h) local k_1, k_2, k_3, k_4; k_1 := `+`(`-`(`*`(`^`(y0, 2)))); k_2 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `*`(k_1)))), 2)))); k_3 := `+`(`-`(`*`(`^`(`+`(y0, `*`(`/`(1, 2), `*`(h, `...
(1)
 

> a := RK(y0, h)
 

`+`(y0, `*`(h, `*`(`+`(`-`(`*`(`/`(1, 6), `*`(`^`(y0, 2)))), `-`(`*`(`/`(1, 3), `*`(`^`(`+`(y0, `-`(`*`(`/`(1, 2), `*`(h, `*`(`^`(y0, 2)))))), 2)))), `-`(`*`(`/`(1, 3), `*`(`^`(`+`(y0, `-`(`*`(`/`(1, ...
`+`(y0, `*`(h, `*`(`+`(`-`(`*`(`/`(1, 6), `*`(`^`(y0, 2)))), `-`(`*`(`/`(1, 3), `*`(`^`(`+`(y0, `-`(`*`(`/`(1, 2), `*`(h, `*`(`^`(y0, 2)))))), 2)))), `-`(`*`(`/`(1, 3), `*`(`^`(`+`(y0, `-`(`*`(`/`(1, ...
(2)
 

Derivatives of the Runge - Kutta operator. 

> simplify(subs([y0 = 1., h = 0.], [seq(diff(a, `$`(h, i)), i = 1 .. 20)]))
 

[-1.000000000, 2.000000000, -6.000000000, 24.00000000, -115.0000000, 600.0000000, -3438.750000, 19530.00000, -107730.0000, 548100.0000, -2442825.000, 9355500.000, -28378350.00, 56756700., -53209406.25...
[-1.000000000, 2.000000000, -6.000000000, 24.00000000, -115.0000000, 600.0000000, -3438.750000, 19530.00000, -107730.0000, 548100.0000, -2442825.000, 9355500.000, -28378350.00, 56756700., -53209406.25...
[-1.000000000, 2.000000000, -6.000000000, 24.00000000, -115.0000000, 600.0000000, -3438.750000, 19530.00000, -107730.0000, 548100.0000, -2442825.000, 9355500.000, -28378350.00, 56756700., -53209406.25...
(3)
 

Derivatives of the actual solution 1/(1+h). 

> simplify(subs([y0 = 1., h = 0.], [seq(diff(`/`(1, `*`(`+`(1, h))), `$`(h, i)), i = 1 .. 20)]))
 

[-1., 2., -6., 24., -120., 720., -5040., 40320., -362880., 3628800., -39916800., 479001600., -6227020800., 0.8717829120e11, -0.1307674368e13, 0.2092278989e14, -0.3556874281e15, 0.6402373706e16, -0.121...
[-1., 2., -6., 24., -120., 720., -5040., 40320., -362880., 3628800., -39916800., 479001600., -6227020800., 0.8717829120e11, -0.1307674368e13, 0.2092278989e14, -0.3556874281e15, 0.6402373706e16, -0.121...
[-1., 2., -6., 24., -120., 720., -5040., 40320., -362880., 3628800., -39916800., 479001600., -6227020800., 0.8717829120e11, -0.1307674368e13, 0.2092278989e14, -0.3556874281e15, 0.6402373706e16, -0.121...
(4)
 

Order 2 operator. 

> RK2 := proc (y0, h) `+`(y0, `-`(`*`(h, `*`(`^`(`+`(y0, `-`(`*`(`/`(1, 2), `*`(h, `*`(`^`(y0, 2)))))), 2))))) end proc; 1
 

proc (y0, h) `+`(y0, `-`(`*`(h, `*`(`^`(`+`(y0, `-`(`*`(`/`(1, 2), `*`(h, `*`(`^`(y0, 2)))))), 2))))) end proc (5)
 

> b := RK2(y0, h); 1
 

`+`(y0, `-`(`*`(h, `*`(`^`(`+`(y0, `-`(`*`(`/`(1, 2), `*`(h, `*`(`^`(y0, 2)))))), 2))))) (6)
 

> simplify(subs([y0 = 1., h = 0.], [seq(diff(b, `$`(h, i)), i = 1 .. 8)]))
 

[-1., 2., -1.500000000, 0, 0, 0, 0, 0] (7)
 

>