381 - A1 - Roundoff

235 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 <module>
    
  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