Runge-Kutta method: Difference between revisions
Content deleted Content added
→{{header|J}}: improve documentation |
Capra Hircus (talk | contribs) m →{{header|Maxima}}: + builtin sol |
||
Line 370: | 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> |
->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}}== |
=={{header|Maxima}}== |
||
<lang maxima> |
<lang maxima>/* Here is how to solve a differential equation */ |
||
'diff(y, x) = x * sqrt(y); |
|||
ode2(%, y, x); |
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], |
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, |
x: x0, |
||
y: y0, |
y: y0, |
||
Line 384: | Line 393: | ||
vy[1]: y0, |
vy[1]: y0, |
||
for i from 1 thru n do ( |
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 |
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])$ |
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> |
for i from 1 step 10 thru 101 do print(x[i], " ", y[i], " ", y[i] - s[i]);</lang> |