Runge-Kutta method: Difference between revisions
Content deleted Content added
m →{{header|C++}}: minor stylistic changes |
m →{{header|REXX}}: added whitespace, optimized a function, simplified some code. |
||
| Line 2,751: | Line 2,751: | ||
╚═══════════════╝ y'(t)=t² √ y(t) ╚═══════════════════════════════════════════╝ |
╚═══════════════╝ y'(t)=t² √ y(t) ╚═══════════════════════════════════════════╝ |
||
</pre> |
</pre> |
||
<lang rexx>/*REXX program uses the Runge─Kutta method to solve the equation: |
<lang rexx>/*REXX program uses the Runge─Kutta method to solve the equation: y'(t) = t² √[y(t)] */ |
||
numeric digits 40; f=digits() % 4 |
numeric digits 40; f= digits() % 4 /*use 40 decimal digs, but only show 10*/ |
||
x0=0; |
x0= 0; x1= 10; dx= .1 /*define variables: X0 X1 DX */ |
||
n=1 + (x1-x0) / dx |
n=1 + (x1-x0) / dx |
||
y.=1; do m=1 for n-1; p=m-1; y.m=RK4(dx, x0 + dx*p, y.p) |
y.=1; do m=1 for n-1; p= m - 1; y.m= RK4(dx, x0 + dx*p, y.p) |
||
end /*m*/ /* [↑] use 4th order Runge─Kutta. */ |
end /*m*/ /* [↑] use 4th order Runge─Kutta. */ |
||
w=digits() % 2 |
w= digits() % 2 /*W: width used for displaying numbers.*/ |
||
say center('X', f, "═") center('Y', w+2, "═") center("relative error", w+8, '═') /*hdr*/ |
say center('X', f, "═") center('Y', w+2, "═") center("relative error", w+8, '═') /*hdr*/ |
||
do i=0 to n-1 by 10; x=(x0 + dx*i) / 1; |
do i=0 to n-1 by 10; x= (x0 + dx*i) / 1; $= y.i / (x*x/4+1)**2 - 1 |
||
say center(x, f) fmt(y.i) left('', 2 + ($>=0) ) fmt($) |
say center(x, f) fmt(y.i) left('', 2 + ($>=0) ) fmt($) |
||
end /*i*/ /*└┴┴┴───◄─────── aligns positive #'s. */ |
end /*i*/ /*└┴┴┴───◄─────── aligns positive #'s. */ |
||
exit /*stick a fork in it, we're all done. */ |
exit /*stick a fork in it, we're all done. */ |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
fmt: z=right( format( |
fmt: parse arg z; z= right( format(z, w, f), w); hasE= pos('E', z)>0; has.= pos(., z)>0 |
||
jus=has. & \hasE; |
jus= has. & \hasE; T= 'T'; if jus then z= left( strip( strip(z, T, 0), T, .), w) |
||
return translate(right(z, (z>=0) + w + 5*hasE + 2*(jus & (z<0) ) ), |
return translate( right(z, (z>=0) + w + 5*hasE + 2*(jus & (z<0) ) ), 'e', "E") |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
RK4: procedure; parse arg dx,x,y; dxH=dx/2; k1= dx * (x ) * sqrt(y ) |
RK4: procedure; parse arg dx,x,y; dxH= dx/2; k1= dx * (x ) * sqrt(y ) |
||
k2= dx * (x + dxH) * sqrt(y + k1/2) |
k2= dx * (x + dxH) * sqrt(y + k1/2) |
||
k3= dx * (x + dxH) * sqrt(y + k2/2) |
k3= dx * (x + dxH) * sqrt(y + k2/2) |
||
k4= dx * (x + dx ) * sqrt(y + k3 ) |
k4= dx * (x + dx ) * sqrt(y + k3 ) |
||
return y + (k1 + k2*2 + k3*2 + k4) / 6 |
return y + (k1 + k2*2 + k3*2 + k4) / 6 |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||