Engineering Computation in Dr. Kyle Horne Department of Mechanical - PowerPoint PPT Presentation

Engineering Computation in Dr. Kyle Horne Department of Mechanical Engineering Spring, 2018

What is It? Explicit Implicit 140 120 100 Position y [px] 15 80 10 60 Position y [px] 5 40 0 20 5 20 40 60 80 100 120 140 240 r =18.80 [µm] Position x [px] Probe 10 235 Laser 15 230 225 15 10 5 0 5 10 15 Position y [µm] Position x [px] 240 r =32.30 [µm] 235 Fit 230 24.5 Temperature T [C] 225 Data 24.0 240 r =44.80 [µm] 23.5 235 230 23.0 225 22.5 230 220 210 200 190 180 170 Position x [µm] 20 25 30 35 40 45 Radial Position r [ um ] Engineering Computation in Python; Dr. Kyle Horne 2 / 44

Where to Learn More? Python Language Numpy (Array Support) Scipy (Numerical Methods) Sympy (Symbolic Math) Matplotlib (Plotting) Engineering Computation in Python; Dr. Kyle Horne 3 / 44

Basic Usage Procedural Execution, Comments and Text Formatting g = 9.81 h = 100.0 # h = 0.5*g*t**2 -> t = sqrt(2*h/g) import math t = math . sqrt (2* h / g ) print( 't = %f [s]' % t ) Results t = 4.515236 Engineering Computation in Python; Dr. Kyle Horne 4 / 44

Basic Usage Lists and Arrays import pylab as pl A = [1.0,2.0,3.0] A . append (4.0) B = pl . array ( A ) C = 2* B print( A ) print( B ) print( C ) Results [1.0 , 2.0 , 3.0 , 4.0] [ 1. 2. 3. 4.] [ 2. 4. 6. 8.] Engineering Computation in Python; Dr. Kyle Horne 5 / 44

Flow Control If-Then-Else a = 1 if a ==1: print( 'a=1' ) elif a ==2: print( 'a=2' ) else: print( 'a!=1' ) Results a=1 Engineering Computation in Python; Dr. Kyle Horne 6 / 44

Flow Control For-Loops s = 0 for k in range (4): x = 2.0* k -1.0 s = s + x print( 'k = %5d; x = %10f' %( k , x )) print( 's = %f' % s ) Results k = 0; x = -1.000000 k = 1; x = 1.000000 k = 2; x = 3.000000 k = 3; x = 5.000000 s = 8.000000 Engineering Computation in Python; Dr. Kyle Horne 7 / 44

Flow Control Functions def f ( x ): y = x **2-1.0 return y g = lambda x : x **2+1.0 x = 2.0 y = f ( x ) print( 'x = %f; y = %f; g(x) = %f' %( x , y , g ( x ))) Results x = 2.000000; y = 3.000000; g(x) = 5.000000 Engineering Computation in Python; Dr. Kyle Horne 8 / 44

Plotting Line Plots and Boilerplate 0.0 import pylab as pl N = 100 0.2 x = pl . linspace (0,1, N ) Ordinate y [-] 0.4 f = lambda x : x **2-1 y = f ( x ) 0.6 fig = pl . figure ( 0.8 tight_layout = True , figsize =(5,4)) 1.0 0.0 0.2 0.4 0.6 0.8 1.0 ax = fig . add_subplot (1,1,1) Abscissa x [-] ax . plot ( x , y , 'b-' , lw =2) ax . plot ( x [::10], y [::10], 'rs' , mew =0) ax . set_xlabel ( 'Abscissa $x$ [-]' ) ax . set_ylabel ( 'Ordinate $y$ [-]' ) pl . savefig ( 'examplePlot.pdf' ) Engineering Computation in Python; Dr. Kyle Horne 9 / 44

Calculation of π (Series) Theory Leibniz’s series formula: � � ∞ ( − 1 ) k ∑ π = 4 2 k + 1 k = 0 3.121595; N = 50 3.8 10 0 Term Magnitude v k v k Terms Sum 3.4 3.0 10 1 2.6 0 10 20 30 40 50 0 10 20 30 40 50 Term Number k Term Count N Engineering Computation in Python; Dr. Kyle Horne 10 / 44

Calculation of π (Series) I Example Solution import numpy 1 import pylab 2 3 N = 50 4 5 terms = numpy . empty ( N ) 6 for k in range ( N ): 7 terms [ k ] = 4.0*(-1)** k /(2.0* k +1.0) 8 9 sums = numpy . empty ( N ) 10 for k in range ( N ): 11 sums [ k ] = terms [: k +1]. sum () 12 13 numbers = numpy . linspace (0, N -1, N ) 14 15 fig = pylab . figure ( tight_layout = True , figsize =(5,4)) 16 ax = fig . add_subplot (1,1,1) 17 ax . plot ( numbers , abs ( terms ), '-' , lw =2) 18 ax . set_yscale ( 'log' ) 19 ax . set_xlabel ( 'Term Number $k$' ) 20 Engineering Computation in Python; Dr. Kyle Horne 11 / 44

Calculation of π (Series) II Example Solution ax . set_ylabel ( 'Term Magnitude $v_k$' ) 21 ax . set_title ( r' ' ) 22 fig . savefig ( 'piTermsPlot.pdf' ) 23 24 fig = pylab . figure ( tight_layout = True , figsize =(5,4)) 25 ax = fig . add_subplot (1,1,1) 26 ax . axhline ( numpy . pi , color = 'k' , linestyle = '--' , linewidth =2) 27 ax . plot ( numbers , sums , '.' , lw =2) 28 ax . set_yticks ([2.6,3.0, numpy . pi ,3.4,3.8]) 29 ax . set_yticklabels ([ '2.6' , '3.0' , '$\pi$' , '3.4' , '3.8' ]) 30 ax . set_xlabel ( r'Term Count $N$' ) 31 ax . set_ylabel ( r'Terms Sum $\sum v_k$' ) 32 ax . text (5,3.8, r'$\pi \approx %f;\,N = %d$' %( terms . sum (), N )) 33 fig . savefig ( 'piSumsPlot.pdf' ) 34 Engineering Computation in Python; Dr. Kyle Horne 12 / 44

Calculation of π (Monte Carlo) 3.07200; N = 500 1.0 A c = π R 2 0.8 A s = ( 2 R ) 2 y -Coordinate 0.6 A c = π A s 4 0.4 � x 2 i + y 2 r i = i 0.2 i ∈ [ 1 , N ] ( r i < 1 ) count   π 0.0 4 = lim   0.0 0.2 0.4 0.6 0.8 1.0 N N → ∞ x -Coordinate Engineering Computation in Python; Dr. Kyle Horne 13 / 44

Calculation of π (Monte Carlo) I Example Solution import pylab 1 import numpy 2 from numpy.random import rand 3 4 def calculatePI ( N ): 5 x = rand ( N ,2) 6 r = numpy . sqrt ( ( x **2). sum (1) ) 7 c = ( r <1.0). sum () 8 9 pi = 4.0* c / float ( N ) 10 11 th = numpy . linspace (0, numpy . pi /2,100) 12 xt = numpy . cos ( th ) 13 yt = numpy . sin ( th ) 14 15 xp = x [:,0] 16 yp = x [:,1] 17 18 fig = pylab . figure ( tight_layout = True , figsize =(4,4)) 19 = fig . add_subplot (1,1,1, aspect =1.0) ax 20 Engineering Computation in Python; Dr. Kyle Horne 14 / 44

Calculation of π (Monte Carlo) II Example Solution ax . plot ( xt , yt , 'k--' , lw =2) 21 ax . plot ( xp [ r <1], yp [ r <1], 'x' , ms =7, mew =1) 22 ax . plot ( xp [ r >1], yp [ r >1], '+' , ms =7, mew =1) 23 ax . set_xlabel ( '$x$-Coordinate' ) 24 ax . set_ylabel ( '$y$-Coordinate' ) 25 ax . set_title ( r'$\pi \approx %.5f; N=%d$' %( pi , N ) ) 26 fig . savefig ( 'randomPI.pdf' ) 27 28 calculatePI (500) 29 Engineering Computation in Python; Dr. Kyle Horne 15 / 44

Sieve of Eratosthenes Theory Sift the Twos and sift the Threes, The Sieve of Eratosthenes! When the multiples sublime, The numbers that remain are Prime. Algorithm ( n ) 10 3 1 P = { 1 , 2 , . . . , N } f ( n ) = n / log ( n ) 1.3 f ( n )/ ( n ) Number of Primes [#] 2 For i ∈ [ 2 , N ] 1.2 10 2 Ratio [1] 1 if P [ i ] is 0, then skip i 1.1 2 For j ∈ [ 2 i , N ] 1.0 10 1 P [ j ] = 0 1 0.9 0.8 10 1 10 2 10 3 10 4 Number n [#] Engineering Computation in Python; Dr. Kyle Horne 16 / 44

Sieve of Eratosthenes I Example Solution import numpy 1 import pylab 2 3 N = 10000+1 4 P = numpy . empty ( N , dtype = int ) 5 6 for i in range ( N ): 7 P [ i ] = i 8 for i in range (2, N ): 9 if P [ i ]==0: 10 continue 11 for j in range (2* i , N , P [ i ]): 12 P [ j ] = 0 13 14 s = 2 15 px = ( P [ P !=0])[ s :] 16 py = numpy . linspace ( s , px . size , px . size ) 17 fy = px / numpy . log ( px ) 18 19 fig = pylab . figure ( tight_layout = True , figsize =(5,4)) 20 Engineering Computation in Python; Dr. Kyle Horne 17 / 44

Sieve of Eratosthenes II Example Solution ax = fig . add_subplot (1,1,1) 21 ax . set_xscale ( 'log' ) 22 ax . set_yscale ( 'log' ) 23 lp , = ax . plot ( px , py , '.' ) 24 lf , = ax . plot ( px , fy , '-' , lw =2) 25 ax . set_xlabel ( 'Number $n$ [#]' ) 26 ax . set_ylabel ( 'Number of Primes $\pi$ [#]' ) 27 ax2 = ax . twinx () 28 lr , = ax2 . plot ( px , fy / py , 'C2_' , mew =2) 29 ax2 . set_ylabel ( r'Ratio [1]' ) 30 lines = [ lp , lf , lr ] 31 labels = [ r'$\pi(n)$' , r'$f(n)=n/log(n)$' , r'$f(n)/\pi(n)$' ] 32 ax2 . legend ( lines , labels , loc = 'upper center' ) 33 fig . savefig ( 'seivePlot.pdf' ) 34 Engineering Computation in Python; Dr. Kyle Horne 18 / 44

Orbital Simulation Theory 1 au Engineering Computation in Python; Dr. Kyle Horne 19 / 44

Orbital Simulation Solver Compute the acceleration on the Earth using Newton’s law of gravitation. Use this to update velocity at each time step. Use velocity to update position at each time step. � t � x ( t ) = � x 0 + � v ( τ ) d τ r 0 ˆ � x 0 = i 0 � t � G · M � � � ˆ v ( t ) = v 0 + a ( τ ) d τ � v 0 = j r 0 0 − G · M r 0 = 1 au � a ( � x 3 � x ) = x The integrals can be converted into much simpler expressions by using Euler’s explicit method: � x + ∆ x f ′ ( ξ ) d ξ ⋍ f ( x ) + ∆ x · f ′ ( x ) f ( x + ∆ x ) = f ( x ) + x Engineering Computation in Python; Dr. Kyle Horne 20 / 44

Orbital Simulation Algorithm Orbit Sun Earth x ∈ R 2 = ( r 0 ) ˆ 1 � i 1.00 �� � 0.75 v ∈ R 2 = ˆ G · M 2 � j r 0 0.50 3 t ∈ R = 0 Position y [au] 0.25 4 While t < t end : 0.00 1 � x = � x + ∆ t · � v 0.25 2 r = | � x | 0.50 a = − G · M 3 � r 3 � x 0.75 4 � v = � v + ∆ t · � a 5 t = t + ∆ t 1.00 1.0 0.5 0.0 0.5 1.0 Position x [au] Engineering Computation in Python; Dr. Kyle Horne 21 / 44

Orbital Simulation I Example Solution import pylab 1 import numpy 2 3 MS = 1.98892e30 4 au = 149.598e9 5 = 6.67E-11 G 6 7 = 365 N 8 9 U0 = numpy . sqrt ( G * MS / au ) 10 Y0 = 31556926.0 11 dt = Y0 /365.0 12 13 xh = numpy . empty ( ( N ,2) ) 14 15 x = numpy . array ([ au ,0.0]) 16 v = numpy . array ([0.0, U0 ]) 17 18 xh [0,:] = x 19 20 Engineering Computation in Python; Dr. Kyle Horne 22 / 44