#Crib: http://www.math.uni-wuppertal.de/~xsc/literatur/PXSCENGL.pdf # INT n = 2; MODE R = LONG LONG REAL; MODE RF = PROC(R)R; RF rsin = long long sin, rcos = long long cos; MODE SV = [n]R; PROC hlf = (SV x)SV: (SV a; FOR j TO n DO a[j] := 0.5*x[j] OD; a); PROC prs = (STRING m, SV x)VOID: (print((m, newline)); FOR j TO n DO print((x[j], newline)) OD); OP * = (R s, SV v)SV: (SV a; FOR j TO n DO a[j] := s*v[j] OD; a); OP * = (REAL s, SV v)SV: (SV a; FOR j TO n DO a[j] := s*v[j] OD; a); OP + = (SV x, y)SV: (SV a; FOR j TO n DO a[j] := x[j] + y[j] OD; a); # This code is to solve the ordinary differential equation dy/dx = f(x, y) # # in particular a circular harmonic oscillator # PROC f = (R tau, SV y)SV: (-y[2], y[1]); PROC step = (R dx, x, SV y)SV: ( SV k1 = dx*f(x,y); SV k2 = dx*f(x+dx/2, y + hlf(k1)); SV k3 = dx*f(x+dx/2, y + hlf(k2)); SV k4 = dx*f(x+dx, y+k3); y + 1.0/R(6)*(k1+2.0*k2 + 2.0*k3 + k4)); SV x := (1, 0); (INT c = 10000; prs("ini", x); TO c DO x := step (1/c, 0, x) OD; print((c, newline)); prs("end", x); print((1.0 - (x[1]*x[1] + x[2]*x[2]), newline, "cos ", x[1]-rcos(1), newline, "sin ", x[2]-rsin(1)))) # ini +1.0000000000000000000000000000000000000000000000000000000e +0 +0.0000000000000000000000000000000000000000000000000000000e +0 +10000 end +5.4030230586813971810212490461842959905274347886258519194e -1 +8.4147098480789650620219196619657541787481134949039919280e -1 +1.3888888871527777777776813368057966579928600000000000000e -22 cos +7.0118829717545299532043305824466296427258410000000000000e -19 sin -4.5031035543372358174775171130797187286931686700000000000e -19 #