Python and Viv · Numerics
Euler against the exact solution
What to notice
errors[-1]counts from the end in both; inside the loop the Viv side readserrors[k].errors.append(e)iserrors = append(errors, e).N0 * 2 ** kisN0 * 2 ^ k, and the step count is a double in Viv; the:4dformat prints it as a whole number.- The parameter
fofeulerand the globaldef fhave the same name in both programs; insideeulerthe parameter is the one called. math.expis the builtinexp. Viv computes it with its own libm, the same in every engine; the printed digits agree with the platform's here.
Python euler_vs_exact.py · 36 lines
Viv euler_vs_exact.viv · 45 lines
py2viv translated from euler_vs_exact.py
What py2viv adds before the program
(nothing)
(nothing)
52 lines: a header and the guard functions its translation calls
1###
2Translated from euler_vs_exact.py by py2viv 0.2.0 (numeric Python -> Viv).
3Python source sha256: 9e3785fe06123be48cce68d41e3e03600cef9d8e4226e7e99931fd530dc287de
4Each print(...) is one print(...) line; fmt(x) writes a float as Python's repr does.
5A `raise`, a failed assert and every guard `stop` with Python's exception line.
6###
7
8#| An int result, as Python's int: a -0.0 becomes 0 (an int has no sign on zero), and the run
9#| stops at 2^53, past which a double no longer holds every integer (Python's int does).
10def py_int(x) {
11 if abs(x) >= 9007199254740992 {
12 stop("py2viv: an integer reached 2^53; Python's int is exact there and a double is not")
13 }
14 return x + 0.0
15}
16
17#| a / b with a float: ZeroDivisionError on a zero divisor.
18def py_div_f(a, b) {
19 if b == 0 {
20 stop("ZeroDivisionError: float division by zero")
21 }
22 return a / b
23}
24
25#| a ** b of two ints, b >= 0: exact (by squaring, every step checked below 2^53).
26def py_ipow(a, b) {
27 if b < 0 {
28 stop("py2viv: an int to a negative power is a float in Python; this translation typed it as an int")
29 }
30 r = 1
31 base = a
32 e = b
33 while e > 0 {
34 if e % 2 == 1 {
35 r = py_int(r * base)
36 }
37 e = floor(e / 2)
38 if e > 0 {
39 base = py_int(base * base)
40 }
41 }
42 return r
43}
44
45#| math.exp(x): OverflowError when a finite x overflows.
46def py_mexp(x) {
47 r = exp(x)
48 if r == 1.0 / 0.0 and x < 1.0 / 0.0 {
49 stop("OverflowError: math range error")
50 }
51 return r
52}
1 imports
4import math
2# nothing to import: exp is a builtin
(nothing)
2 the problem: y' = -2 t y, y(0) = 1, exact y = exp(-t^2)
7def f(t, y):
8 return -2.0 * t * y
9
10def exact(t):
11 return math.exp(-t * t)
5def f(t, y) {
6 return -2.0 * t * y
7}
8
9def exact(t) {
10 return exp(-t * t)
11}
58def f(t, y) {
59 return -2.0 * t * y
60}
61def exact(t) {
62 return py_mexp(-t * t)
63}
3 forward Euler
14def euler(f, y0, t_end, n):
15 h = t_end / n
16 t, y = 0.0, y0
17 for _ in range(n):
18 y += h * f(t, y)
19 t += h
20 return y
14def euler(f, y0, t_end, n) {
15 h = t_end / n
16 t, y = 0.0, y0
17 for k to n {
18 y += h * f(t, y)
19 t += h
20 }
21 return y
22}
65def euler(f, y0, t_end, n) {
66 h = py_div_f(t_end, n)
67 t = 0.0
68 y = y0
69 for _ to n {
70 y += h * f(t, y)
71 t += h
72 }
73 return y
74}
4 the error halves when the step halves
23T_END = 2.0
24N0 = 10
25REFINEMENTS = 5
26errors = []
27for k in range(REFINEMENTS):
28 n = N0 * 2 ** k
29 y = euler(f, 1.0, T_END, n)
30 errors.append(abs(y - exact(T_END)))
31 print(f"{n:4d} steps y(2) = {y:.9f} error = {errors[-1]:.3e}")
32for k in range(1, REFINEMENTS):
33 print(f"halving the step divides the error by {errors[k - 1] / errors[k]:.3f}")
25T_END = 2.0
26N0 = 10
27REFINEMENTS = 5
28errors = []
29for k to REFINEMENTS {
30 n = N0 * 2 ^ k
31 y = euler(f, 1.0, T_END, n)
32 errors = append(errors, abs(y - exact(T_END)))
33 print(f"{n:4d} steps y(2) = {y:.9f} error = {errors[k]:.3e}")
34}
35for k from 1 to REFINEMENTS {
36 print(f"halving the step divides the error by {errors[k - 1] / errors[k]:.3f}")
37}
76T_END = 2.0
77N0 = 10
78REFINEMENTS = 5
79errors = []
80for k to REFINEMENTS {
81 n = py_int(N0 * py_ipow(2, k))
82 y = euler(f, 1.0, T_END, n)
83 errors = append(errors, abs(y - exact(T_END)))
84 print(fmt(n + 0.0, "4d"), "steps y(2) =", fmt(y, ".9f"), " error =", fmt(errors[-1], ".3e"))
85}
86for k from 1 to REFINEMENTS {
87 print("halving the step divides the error by", fmt(py_div_f(errors[py_int(k - 1)], errors[k]), ".3f"))
88}
5 checks
36assert abs(errors[-2] / errors[-1] - 2.0) < 0.1
40test "first order" {
41 assert_near(errors[-2] / errors[-1], 2.0, 0.1)
42}
90if not (abs(py_div_f(errors[-2], errors[-1]) - 2.0) < 0.1) {
91 stop("AssertionError")
92}
What they print
identical — the two programs print the same text.
Dimmed lines are left out of the comparison. Each Viv test block prints one PASS line; the Python checks the same facts with assert, which prints nothing.
$ python3 euler_vs_exact.py
10 steps y(2) = 0.005526592 error = 1.279e-0220 steps y(2) = 0.012023052 error = 6.293e-0340 steps y(2) = 0.015213620 error = 3.102e-0380 steps y(2) = 0.016776771 error = 1.539e-03160 steps y(2) = 0.017549325 error = 7.663e-04halving the step divides the error by 2.032halving the step divides the error by 2.029halving the step divides the error by 2.016halving the step divides the error by 2.008
$ viv run euler_vs_exact.viv
10 steps y(2) = 0.005526592 error = 1.279e-0220 steps y(2) = 0.012023052 error = 6.293e-0340 steps y(2) = 0.015213620 error = 3.102e-0380 steps y(2) = 0.016776771 error = 1.539e-03160 steps y(2) = 0.017549325 error = 7.663e-04halving the step divides the error by 2.032halving the step divides the error by 2.029halving the step divides the error by 2.016halving the step divides the error by 2.008PASS first order
$ viv run euler_vs_exact.py (py2viv’s translation)
The translation prints exactly what the Python prints.
10 steps y(2) = 0.005526592 error = 1.279e-0220 steps y(2) = 0.012023052 error = 6.293e-0340 steps y(2) = 0.015213620 error = 3.102e-0380 steps y(2) = 0.016776771 error = 1.539e-03160 steps y(2) = 0.017549325 error = 7.663e-04halving the step divides the error by 2.032halving the step divides the error by 2.029halving the step divides the error by 2.016halving the step divides the error by 2.008