Runge-Kutta method: Difference between revisions
Content added Content deleted
(→{{header|Python}}: removing solution involving nested lambdas) |
(Added solution for Action!) |
||
Line 65: | Line 65: | ||
9.0 451.56246 -0.00004072 |
9.0 451.56246 -0.00004072 |
||
10.0 675.99995 -0.00005098 |
10.0 675.99995 -0.00005098 |
||
</pre> |
|||
=={{header|Action!}}== |
|||
Calculations on a real Atari 8-bit computer take quite long time. It is recommended to use an emulator capable with increasing speed of Atari CPU. |
|||
{{libheader|Action! Tool Kit}} |
|||
{{libheader|Action! Real Math}} |
|||
<lang Action!>INCLUDE "D2:PRINTF.ACT" ;from the Action! Tool Kit |
|||
INCLUDE "H6:REALMATH.ACT" |
|||
DEFINE PTR="CARD" |
|||
REAL one,two,four,six |
|||
PROC Init() |
|||
IntToReal(1,one) |
|||
IntToReal(2,two) |
|||
IntToReal(4,four) |
|||
IntToReal(6,six) |
|||
RETURN |
|||
PROC Fun=*(REAL POINTER x,y,res) |
|||
DEFINE JSR="$20" |
|||
DEFINE RTS="$60" |
|||
[JSR $00 $00 ;JSR to address set by SetFun |
|||
RTS] |
|||
PROC SetFun(PTR p) |
|||
PTR addr |
|||
addr=Fun+1 ;location of address of JSR |
|||
PokeC(addr,p) |
|||
RETURN |
|||
PROC Rate(REAL POINTER x,y,res) |
|||
REAL tmp |
|||
Sqrt(y,tmp) ;tmp=sqrt(y) |
|||
RealMult(x,tmp,res) ;res=x*sqrt(y) |
|||
RETURN |
|||
PROC RK4(PTR f REAL POINTER dx,x,y,res) |
|||
REAL k1,k2,k3,k4,dx2,k12,k22,tmp1,tmp2,tmp3 |
|||
SetFun(f) |
|||
Fun(x,y,tmp1) ;tmp1=f(x,y) |
|||
RealMult(dx,tmp1,k1) ;k1=dx*f(x,y) |
|||
RealDiv(dx,two,dx2) ;dx2=dx/2 |
|||
RealDiv(k1,two,k12) ;k12=k1/2 |
|||
RealAdd(x,dx2,tmp1) ;tmp1=x+dx/2 |
|||
RealAdd(y,k12,tmp2) ;tmp2=y+k1/2 |
|||
Fun(tmp1,tmp2,tmp3) ;tmp3=f(x+dx/2,y+k1/2) |
|||
RealMult(dx,tmp3,k2) ;k2=dx*f(x+dx/2,y+k1/2) |
|||
RealDiv(k2,two,k22) ;k22=k2/2 |
|||
RealAdd(y,k22,tmp2) ;tmp2=y+k2/2 |
|||
Fun(tmp1,tmp2,tmp3) ;tmp3=f(x+dx/2,y+k2/2) |
|||
RealMult(dx,tmp3,k3) ;k3=dx*f(x+dx/2,y+k2/2) |
|||
RealAdd(x,dx,tmp1) ;tmp1=x+dx |
|||
RealAdd(y,k3,tmp2) ;tmp2=y+k3 |
|||
Fun(tmp1,tmp2,tmp3) ;tmp3=f(x+dx,y+k3) |
|||
RealMult(dx,tmp3,k4) ;k4=dx*f(x+dx,y+k3) |
|||
RealAdd(k2,k3,tmp1) ;tmp1=k2+k3 |
|||
RealMult(two,tmp1,tmp2) ;tmp2=2*k2+2*k3 |
|||
RealAdd(k1,tmp2,tmp1) ;tmp3=k1+2*k2+2*k3 |
|||
RealAdd(tmp1,k4,tmp2) ;tmp2=k1+2*k2+2*k3+k4 |
|||
RealDiv(tmp2,six,tmp1) ;tmp1=(k1+2*k2+2*k3+k4)/6 |
|||
RealAdd(y,tmp1,res) ;res=y+(k1+2*k2+2*k3+k4)/6 |
|||
RETURN |
|||
PROC Calc(REAL POINTER x,res) |
|||
REAL tmp1,tmp2 |
|||
RealMult(x,x,tmp1) ;tmp1=x*x |
|||
RealDiv(tmp1,four,tmp2) ;tmp2=x*x/4 |
|||
RealAdd(tmp2,one,tmp1) ;tmp1=x*x/4+1 |
|||
Power(tmp1,two,res) ;res=(x*x/4+1)^2 |
|||
RETURN |
|||
PROC RelError(REAL POINTER a,b,res) |
|||
REAL tmp |
|||
RealDiv(a,b,tmp) ;tmp=a/b |
|||
RealSub(tmp,one,res) ;res=a/b-1 |
|||
RETURN |
|||
PROC Main() |
|||
REAL x0,x1,x,dx,y,y2,err,tmp1,tmp2 |
|||
CHAR ARRAY s(20) |
|||
INT i,n |
|||
Put(125) PutE() ;clear the screen |
|||
MathInit() |
|||
Init() |
|||
PrintF("%-2S %-11S %-8S%E","x","y","rel err") |
|||
IntToReal(0,x0) |
|||
IntToReal(10,x1) |
|||
ValR("0.1",dx) |
|||
RealSub(x1,x0,tmp1) ;tmp1=x1-x0 |
|||
RealDiv(tmp1,dx,tmp2) ;tmp2=(x1-x0)/dx |
|||
n=RealToInt(tmp2) ;n=(x1-x0)/dx |
|||
i=0 |
|||
IntToReal(1,y) |
|||
DO |
|||
IntToReal(i,tmp1) ;tmp1=i |
|||
RealMult(dx,tmp1,tmp2) ;tmp2=i*dx |
|||
RealAdd(x0,tmp2,x) ;x=x0+i*dx |
|||
IF i MOD 10=0 THEN |
|||
Calc(x,y2) |
|||
RelError(y,y2,err) |
|||
StrR(x,s) PrintF("%-2S ",s) |
|||
StrR(y,s) PrintF("%-11S ",s) |
|||
StrR(err,s) PrintF("%-8S%E",s) |
|||
FI |
|||
i==+1 |
|||
IF i>n THEN EXIT FI |
|||
RK4(rate,dx,x,y,tmp1) ;tmp1=rk4(rate,dx,x0+dx*(i-1),y) |
|||
RealAssign(tmp1,y) ;y=rk4(rate,dx,x0+dx*(i-1),y) |
|||
OD |
|||
RETURN</lang> |
|||
{{out}} |
|||
[https://gitlab.com/amarok8bit/action-rosetta-code/-/raw/master/images/Runge-Kutta_method.png Screenshot from Atari 8-bit computer] |
|||
<pre> |
|||
x y rel err |
|||
0 1 0 |
|||
1 1.56249977 -1.3E-07 |
|||
2 3.99999882 -2.9E-07 |
|||
3 10.56249647 -2.9E-07 |
|||
4 24.99999228 -2.9E-07 |
|||
5 52.56248607 -2.0E-07 |
|||
6 99.99997763 -2.1E-07 |
|||
7 175.562459 -1.8E-07 |
|||
8 288.999935 -1.9E-07 |
|||
9 451.562406 0 |
|||
10 675.999869 -1.4E-07 |
|||
</pre> |
</pre> |
||