40 lines
1 KiB
Text
40 lines
1 KiB
Text
/* Here is how to solve a differential equation */
|
|
'diff(y, x) = x * sqrt(y);
|
|
ode2(%, y, x);
|
|
ic1(%, x = 0, y = 1);
|
|
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)),
|
|
x: x0,
|
|
y: y0,
|
|
vx: makelist(0, n + 1),
|
|
vy: makelist(0, n + 1),
|
|
vx[1]: x0,
|
|
vy[1]: y0,
|
|
for i from 1 thru n do (
|
|
k1: bfloat(h * f(x, y)),
|
|
k2: bfloat(h * f(x + h / 2, y + k1 / 2)),
|
|
k3: bfloat(h * f(x + h / 2, y + k2 / 2)),
|
|
k4: bfloat(h * f(x + h, y + k3)),
|
|
vy[i + 1]: y: y + (k1 + 2 * k2 + 2 * k3 + k4) / 6,
|
|
vx[i + 1]: x: x + h
|
|
),
|
|
[vx, vy]
|
|
)$
|
|
|
|
[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)$
|
|
|
|
for i from 1 step 10 thru 101 do print(x[i], " ", y[i], " ", y[i] - s[i]);
|