Runge-Kutta method: Difference between revisions

Line 370:

->y[{1,2,3,4,5,6,7,8,9,10}]->{25/16, 4, 169/16, 25, 841/16, 100, 2809/16, 289, 7225/16, 676}</lang>

=={{header|Maxima}}==

<lang maxima>~~'diff(y,x)~~ = x*~~sqrt(y);~~

<lang maxima>/* Here is how to solve a differential equation */

'diff(y, x) = x * sqrt(y);

ode2(%, y, x);

ic1(%, x=0, y=1);

ic1(%, x = 0, y = 1);

factor(solve(%, y)); /* [y=(x^2+4)^2/16] */

factor(solve(%, y)); /* [y = (x^2 + 4)^2 / 16] */

/* The Runge-Kutta solver is builtin */

load(dynamics)$

sol: rk(t * sqrt(y), y, 1, [t, 0, 10, 1.0])$

plot2d([discrete, sol])$

/* An implementation of RK4 for one equation */

rk4(f, x0, y0, x1, n) := block([h, x, y, vx, vy, k1, k2, k3, k4],

h: bfloat((x1 - x0)/(n - 1)),

h: bfloat((x1 - x0) / (n - 1)),

x: x0,

y: y0,

Line 384:

Line 393:

vy[1]: y0,

for i from 1 thru n do (

k1: bfloat(h*f(x, y)),

k1: bfloat(h * f(x, y)),

k2: bfloat(h*f(x + h/2, y + k1/2)),

k2: bfloat(h * f(x + h / 2, y + k1 / 2)),

k3: bfloat(h*f(x + h/2, y + k2/2)),

k3: bfloat(h * f(x + h / 2, y + k2 / 2)),

k4: bfloat(h*f(x + h, y + k3)),

k4: bfloat(h * f(x + h, y + k3)),

vy[i + 1]: y: y + (k1 + 2*k2 + 2*k3 + k4)/6,

vy[i + 1]: y: y + (k1 + 2 * k2 + 2 * k3 + k4) / 6,

vx[i + 1]: x: x + h

),

Line 394:

Line 403:

)$

[x, y]: rk4(lambda([x, y], x*sqrt(y)), 0, 1, 10, 101)$

[x, y]: rk4(lambda([x, y], x * sqrt(y)), 0, 1, 10, 101)$

plot2d([discrete, x, y])$

s: map(lambda([x], (x^2 + 4)^2/16), x)$

s: map(lambda([x], (x^2 + 4)^2 / 16), x)$

for i from 1 step 10 thru 101 do print(x[i], " ", y[i], " ", y[i] - s[i]);</lang>