# 381 - A1 - Roundoff

## 373 days ago by Professor381

eps = 0.0001
# Solving easy "stable" difference equation from sympy import Function, rsolve,Poly from sympy.abc import n y=Function('y') f = y(n + 1) - 2/3*y(n) print("General solution is ",rsolve(f, y(n))) inits={y(0):1} soln = rsolve(f, y(n),inits) print(soln) # points([(n,soln(n)) for n in range(1,10)])
 ('General solution is ', (2/3)**n*C0) (2/3)**n ('General solution is ', (2/3)**n*C0) (2/3)**n
# Solving second order "stable" difference equation from sympy import Function, rsolve,Poly from sympy.abc import n y=Function('y') f = y(n + 2) - 7/6*y(n+1) + 1/3*y(n) print("General solution is ",rsolve(f, y(n))) inits={y(0):1,y(1):2/3} print(rsolve(f, y(n),inits)) eps = 0.01 inits={y(0):1,y(1):2/3+eps} rsolve(f, y(n),inits)
 ('General solution is ', (2/3)**n*C1 + 2**(-n)*C0) (2/3)**n  ('General solution is ', (2/3)**n*C1 + 2**(-n)*C0) (2/3)**n 
# Solving second order "relatively stable" difference equation from sympy import Function, rsolve,Poly from sympy.abc import n y=Function('y') f = y(n + 2) - 17/12*y(n+1) + 1/2*y(n) print("General solution is ",rsolve(f, y(n))) inits={y(0):1,y(1):2/3} rsolve(f, y(n),inits) inits={y(0):1,y(1):2/3+eps} rsolve(f, y(n),inits)
 ('General solution is ', (2/3)**n*C0 + (3/4)**n*C1)  ('General solution is ', (2/3)**n*C0 + (3/4)**n*C1) 
# Solving second order "unstable" difference equation from sympy import Function, rsolve,Poly from sympy.abc import n y=Function('y') f = y(n + 2) - 8/3*y(n+1) + 4/3*y(n) print("General solution is ",rsolve(f, y(n))) inits={y(0):1,y(1):2/3} print(rsolve(f, y(n),inits)) inits={y(0):1,y(1):2/3+eps} rsolve(f, y(n),inits)
 ('General solution is ', 2**n*C0 + (2/3)**n*C1) (2/3)**n  ('General solution is ', 2**n*C0 + (2/3)**n*C1) (2/3)**n 

Now, change to octave mode to evaluate the remainder of these done numerically.

for n = 1:100 exact(n) = (2/3)^(n-1); end p(1) = 1; p(2) = 2/3; for n=3:100 p(n) = 2/3*p(n-1); relerr(n) = abs(p(n)-exact(n))/exact(n); end # plot(p) plot(relerr)
 Traceback (click to the left of this block for traceback) ... ValueError: Argument names should be valid python identifiers. Traceback (most recent call last): p(n) = 2/3*p(n-1); File "", line 1, in File "/home/sageserver/sage-8.7/local/lib/python2.7/site-packages/sagenb/misc/support.py", line 685, in preparse_worksheet_cell s = preparse_file(s, globals=globals) File "/home/sageserver/sage-8.7/local/lib/python2.7/site-packages/sage/repl/preparse.py", line 1343, in preparse_file numeric_literals=not numeric_literals)) File "/home/sageserver/sage-8.7/local/lib/python2.7/site-packages/sage/repl/preparse.py", line 1236, in preparse L = preparse_calculus(L) File "/home/sageserver/sage-8.7/local/lib/python2.7/site-packages/sage/repl/preparse.py", line 939, in preparse_calculus raise ValueError("Argument names should be valid python identifiers.") ValueError: Argument names should be valid python identifiers.
for n = 1:100 exact(n) = (2/3)^(n-1); end p(1) = 1; p(2) = 2/3; for n=3:100 p(n) = 7/6*p(n-1)-1/3*p(n-2); relerr(n) = abs(p(n)-exact(n))/exact(n); end plot(p) plot(relerr)
  8e-14 |-------------------------------------------------------------------| | + + + + | 7e-14 |-+ +-| | | | | 6e-14 |-+ ****| | *** | 5e-14 |-+ ***** +-| | *** | | ***** | 4e-14 |-+ ***** +-| | ****** | 3e-14 |-+ ***** +-| | ******* | | ***** | 2e-14 |-+ ***** +-| | **** | 1e-14 |-+ *** +-| | **** | | ***** + + + + | 0 |-------------------------------------------------------------------| 0 20 40 60 80 100   8e-14 |-------------------------------------------------------------------| | + + + + | 7e-14 |-+ +-| | | | | 6e-14 |-+ ****| | *** | 5e-14 |-+ ***** +-| | *** | | ***** | 4e-14 |-+ ***** +-| | ****** | 3e-14 |-+ ***** +-| | ******* | | ***** | 2e-14 |-+ ***** +-| | **** | 1e-14 |-+ *** +-| | **** | | ***** + + + + | 0 |-------------------------------------------------------------------| 0 20 40 60 80 100 
for n = 1:100 exact(n) = (2/3)^(n-1); end p(1) = 1; p(2) = 2/3; for n=3:100 p(n) = 17/12*p(n-1)-1/2*p(n-2); relerr(n) = abs(p(n)-exact(n))/exact(n); end plot(p) plot(relerr)
  1e-09 |-------------------------------------------------------------------| | + + + + | | | | | 8e-10 |-+ +-| | | | | | *| 6e-10 |-+ +*| | *| | * | | * | 4e-10 |-+ *+-| | * | | * | | ** | 2e-10 |-+ ** +-| | ** | | *** | | + + + ****** | 0 |-------------------------------------------------------------------| 0 20 40 60 80 100   1e-09 |-------------------------------------------------------------------| | + + + + | | | | | 8e-10 |-+ +-| | | | | | *| 6e-10 |-+ +*| | *| | * | | * | 4e-10 |-+ *+-| | * | | * | | ** | 2e-10 |-+ ** +-| | ** | | *** | | + + + ****** | 0 |-------------------------------------------------------------------| 0 20 40 60 80 100 
for n = 1:100 exact(n) = (2/3)^(n-1); end p(1) = 1; p(2) = 2/3; for n=3:100 p(n) = 8/3*p(n-1)-4/3*p(n-2); relerr(n) = abs(p(n)-exact(n))/exact(n); end plot(p) plot(relerr)
  3e+29 |-----------------------------------------------------------------| | + + + + | | | 2.5e+29 |-+ +-| | | | | | | 2e+29 |-+ +*| | *| | *| 1.5e+29 |-+ +*| | *| | *| 1e+29 |-+ +*| | *| | *| | *| 5e+28 |-+ +*| | *| | + + + + * | 0 |-----------------------------------------------------------------| 0 20 40 60 80 100   3e+29 |-----------------------------------------------------------------| | + + + + | | | 2.5e+29 |-+ +-| | | | | | | 2e+29 |-+ +*| | *| | *| 1.5e+29 |-+ +*| | *| | *| 1e+29 |-+ +*| | *| | *| | *| 5e+28 |-+ +*| | *| | + + + + * | 0 |-----------------------------------------------------------------| 0 20 40 60 80 100 

Let's put these all together on one graph.

# plotting relative errors for stable, relatively stable, and unstable difference equations top = 25 ns = 1:top; for n = 1:top exact(n) = (2/3)^(n-1); end p1(1) = 1; p1(2) = 2/3; for n=3:top p1(n) = 2/3*p1(n-1); relerr1(n) = abs(p1(n)-exact(n))/exact(n); end p2(1) = 1; p2(2) = 2/3; for n=3:top p2(n) = 7/6*p2(n-1)-1/3*p2(n-2); relerr2(n) = abs(p2(n)-exact(n))/exact(n); end p3(1) = 1; p3(2) = 2/3; for n=3:top p3(n) = 17/12*p3(n-1)-1/2*p3(n-2); relerr3(n) = abs(p3(n)-exact(n))/exact(n); end p4(1) = 1; p4(2) = 2/3; for n=3:top p4(n) = 8/3*p4(n-1)-4/3*p4(n-2); relerr4(n) = abs(p4(n)-exact(n))/exact(n); end plot(ns,relerr1) hold on plot(ns,relerr2) plot(ns,relerr3) plot(ns,relerr4) axis([1, top, 0, 0.0000000000000001])
 top = 25 error: __plt2vv__: vector lengths must match error: called from __plt__>__plt2vv__ at line 482 column 5 __plt__>__plt2__ at line 242 column 14 __plt__ at line 107 column 18 plot at line 223 column 10 /home/sageserver/.sage/temp/travis-PowerEdge-T430/18295/interface/tmp183\ 15 at line 32 column 1 1 |-+ | | | 0.8 |-+ | | | 0.6 |-+ | | | 0.4 |-+ | | | 0.2 |-+ | | | + + + + + 0 |---------------------------------------------------------------------- 0 0.2 0.4 0.6 0.8 1  top = 25 error: __plt2vv__: vector lengths must match error: called from __plt__>__plt2vv__ at line 482 column 5 __plt__>__plt2__ at line 242 column 14 __plt__ at line 107 column 18 plot at line 223 column 10 /home/sageserver/.sage/temp/travis-PowerEdge-T430/18295/interface/tmp18315 at line 32 column 1 1 |-+ | | | 0.8 |-+ | | | 0.6 |-+ | | | 0.4 |-+ | | | 0.2 |-+ | | | + + + + + 0 |---------------------------------------------------------------------- 0 0.2 0.4 0.6 0.8 1 

New experimant with square rooting and then squaring back the same number of times.  You should get what you started with...

x0 = 200 x = x0 n = 50 for k=1:n x = sqrt(x); end for k=1:n x = x^2; end x (x0-x)/x0
 x0 = 200 x = 200 n = 50 x = 190.566 ans = 0.0471687 x0 = 200 x = 200 n = 50 x = 190.566 ans = 0.0471687
# Now do it in Sage mode x0 = 2 x = x0 n = 50 for k in range(n): x = sqrt(x); end for k in range(n): x = x^2; end print(x) print((x0-x)/x0)
 2 0 2 0