Contents: Computer Science and Engineering Stream Lab 1 Lab 2 Lab 3. Lab 4. Lab 5. Lab 6 Lab 7. Lab 8. Lab 9. Lab 10. Programme to compute area, volume and center of gravity. Evaluation of improper integrals , Beta and Gamma functions. Finding gradient, divergent, curl and their geometrical interpretation ‘Computation of basis and dimension for a vector space and graphical representation of linear transformation Computing the inner product and orthogonality Solution of algebraic and transcendental equation by Rogula-Falsi and Newton- Raphson method Interpolation /Extrapolation using Newton’s forward and backward difference for- mula Computation of area under the curve using Trapezoidal, Simpson's (2) and Simp- sons (3)"" rule Solution of ODE of first order and first degree by Taylor's series and Modified Euler’s method Solution of ODE of first order and first degree by Runge-Kutta 4th order method and Milne’s predictor and corrector methodLAB 1: Programme to compute area, volume and cen- ter of gravity 1.1 Objectives: Use python 1. to evaluate double integration. 2. to compute area and volume. 3. to calculate center of gravity of 2D object. Syntax for the commands used: 1. Data pretty printer in Python. pprint () 2. integrate: integrate(function (variable, min limit, max_limit)) 1.2 Double and triple integration Example 1: Le Evaluate the integral [ [(2? + y?)dydx aa from sympy import x.y, z=symbols ("x y 2") wirintegrate(xs#2ty#=2,(y,0,x) ,(x,0,1)) print (wt) 1/3 Example 2: gacetsy Evaluate the integral f ff (xyz)dzdyde aaa from sympy import + ymbo Cx") ymbol('y') zeSymbol('z') w2= integrate ((xry#z) ,(2,0,3-x-y) ,(y,0,3-x) (0,39) print (w2) 81/80Example 3: Prove that f f(x? + y*)dydx = f [(2? + y%)drdy from sympy import + ymbol ('x') ymbol('y') symbol ('z') WSs integrate (x*#2tys=2,y,x) pprint (#3) wa= integrate (x*#2+y*=2,x,y) pprint (w4) 1.3 Area and Volume Area of the region R in the cartesian form is ff drdy R Example 4: a (fave Find the area of an ellipse by double integration. A=4 ff) dydx a from sympy import + ymbol (x!) symbol ('y') #axSymbol('a') #b=Symbol ('b') u3-44 integrate (1, (y,0,(b/a)=sart (are2-x**2)), (x,0,a)) print (w3) 24. 0*pi 1.4. Area of the region R in the polar form is f [rdrd9 R Example 5: Find the area of the eardioid r = a(1 + cos) by double integration fron sympy import + r=Symbol ('r') symbol ('t') a>Symbol('a') wana w3-2+ integrate (r,(r,0,as(1+cos(t))) ,(t,0,pi)) pprint (w3)1.5 Volume of a solid is given by [ f { dedydz i Example 6: Find the volume of the tetrahedron bounded by the planes x: from sympy import + x=Symbol ("x") yoSymbol ('y") z=Symbol('z') a-Symbol('a') b=Symbol('b') Symbol ('c') w2=integrate(1,(z,0,c%(1-x/a-y/b)) ,(y,0,b4(1-x/a)), (x,0,a)) print (w2) atbtc/6 1.6 Center of Gravity Find the conter of gravity of cardioid . Plot the graph of cardioid and mark the center of gravity. import numpy as np import matplotlib.pyplot as pit import math from sympy import * r=Symbol ('r') teSymbo1('t') a=Symbol ('a') Ti=integrate (cos(t)*r4=2,(r,0,a*(1+cos(t))) ,(t,-pi pid) I2=integrate(r,(r,0,a*(1*cos(t))) ,(t,~pi,pi)) Te11/12 print (1) I=1. subs (2,5) plt.axes(projection = ‘polar') ans rad = np-arange(0, (2 * mp-pi), 0.01) # plotting the cardioid for i in rad r= a + (arap.cos (i) pit. polar(i,r,'g.") plt.polar(0,I,'r.') plt.show() 5¥a/61.7 1 270° Exercise: Evaluate j fee + y)dyde Ans: 0.5 og(2) 2 = >toa(y) Find the fff (7 )dzdyde Ans: -0. 2627 Find the area of positive quadrant of the circle 2? + y? = 16 Ans: 4 Find the volume of the tetrahedron bounded by the planes x=! g+¥sae1 Ans: 4 58LAB 2: Evaluation of improper integrals, Beta and Gamma functions 2.1 Objectives: Use python 1. to evaluate improper integrals using Beta function. 2. to evaluate improper integrals using Gamma function Syntax for the commands used: 1. gamma math. gamma (x) Parameters : x : The number whose gamma value needs to be computed. 2. beta math. beta(x,y) Parameters : x ,y: The numbers whose beta value needs to be computed. 3. Note: We can evaluate improper integral involving infinity by using inf Example 1: Evaluate f e*de. from sympy import * ymbols ('x') wisintegrate (exp (=x) ,(x,0, float (‘inf '))) print (simplify (w1)) Gamma funetion is x(n) = [2° e*2""'de Example 2: Evaluate (5) by using definition from sympy import + yabols('x') wizintegrate (exp (-x)*xe4,(x,0, float ('inf'))) print (simplify (w1)) 24Example 3: Evaluate f e~*‘cos(4t)dt . That is Laplace transform of cos( 4) a fron sympy izport + t,sesymbols('t,s') #' for infinity in sympy ve use oo wis integrate (exp(-s+t)+cos (4t) ,(€,0,00)) display (simplify (w1)) 2a16 [cos (i) dt otherwise a for 2arg (s)| R? geometrically. Find the image of vector (10,0) when it is reflected about y axis import mumpy as up import matplotlib.pyplot as plt v= np.array({{10,0]]) origin = np.array({[0, 0, 0],(0, 0, 0]]) # origin point Aenp.matrix({(-1,0],(0,1]]) Vienp-matrix(V) V2=Anp. transpose (V1) V2enp.array (V2) plt.quiver(+origin, V[:,0], V[:,1], color=['b'], seale=50) plt.quiver(*origin, V2(0,-), V2(1,:], color=['r'], scale=80) plt. show () 0.08 002 0.00 ———$ 002 0.08 “00s 0020002 Another example. B= np.array({[-1,0],(0,1]1) B_coords = B@coords X.LT2 = Bcoords 0, :] y.LT2 = Bleoords{1,:] # Create the figure and axes objects fig, ax = plt.subplots() # Plot the points. x and y are original vectors, x LT1 and y_LT1 are images ax.plot (x,y, 'r0") ax. plot (xLT2,y.LT2,'bo") # Connect the points by lines ax.plot(x,y,'r!,ls="=-") ax. plot (x_LT2,y_LT2,"") # Edit some settings |ax-axvline (x=0, color " 67ax. axhline(y=0,color="k",1s= ax. grid(True) ax.axis((-2,2,-1,21) ax.set_aspect (‘equal’) ax.set_title ("Reflection"); ” Reflection 2.0 : 15 10 os “15-10 -05 00 O05 10 15 20 4.4.3 Rotation: Represent the rotation transformation 7: R? + R? geometrically. Find the image of vector (10,0) when it is rotated by 7/2 radians. import numpy as np import matplotlib.pyplot as plt V = np.array([[10,01]) origin = np.array({[0, 0, 0),(0, 0, 0)]) # origin point A=np.matrix({(0,-1],{1,1]]) Vienp-matrix(V) V2=A*np. transpose (V1) v2=np. array (V2) plt.quiver(+origin, V[:,0], V[:,1], color=['b'], scale=50) plt.quiver(+origin, V2(0,-), V2it,:], color=['r'], scale=S0) plt.show() 68008 02 0.00 0.02 0.08 moos 9.02 9999908 Another example. theta = pi/é R = np.array({[cos (theta) ,-sin(theta)] , (sin(theta),cos(theta)])) R.coords = Récoords XLT = Rcoords (0, ylLT3 = R_coords[t, # Create the figure and axes objects fig, ax = plt.subplots() # Plot the points. x and y are original vectors, x LT1 and y_LT1 are images ax.plot (x,y, !r0") ax. plot (x_LT3,y_LT3,"bo') # Connect the points by lines ax.plot(x,y,"r",1s="==") ax.plot (x_LT3,y/LT3,"b! # Edit some settings ax.axvline(x=0,color="k",1s=":") ax. axhline (y=0,color="k",1s=":" ax. grid (True) ax.axis([-2,2,-1,21) ax.set_aspect (‘equal’) 692.0 1s oS 10 05 05 10 15 20 4.4.4 Shear ‘Transformation Represent the Shear transformation 7: R? + R? geometrically. Find the image of (2,3) under shear transformation. import numpy as np import matplotiib.pyplot as pit V = np-array(([2,3]]) origin = np.array(((0, 0, 0),(0, Aenp.matrix({(1,21,[0,1]]) Vienp. matrix (V) V2=A+np. transpose (V1) Va-np. array (V2) 0))) # origin point print("Image of given vectors is:", V2) plt-quiver(+origin, V[:,0], V{:,1], color=['b'], plt.quiver(+origin, V2(0,:], V2(t,:], color=('r'], scale=20) pit. show ()008 02 0.00 0.02 0.08 “00s 00200 Another example. S = np.array(((1,2],[0,1]]) S_coords = S@coords x_LT4 y_LT4 S_coords[0,:) S_coords[1,:) # Create the figure and axes objects £ig, ax = plt. subplots () # Plot the points. x and y are original vectors, x_LT1 and y_LT1 are images ax.plot (x,y, 'ro') ax. plot (x_LT4,y_LT4, bo") # Connect the points by lines ax.plot (x,y, 'r' yl ax. plot (x_LT4,y_LT4,"b") # Edit some settings ax. axvline (x=0,color="k",1s=":") ax. axhline (y=0, color ax. grid (True) ax.axis((-2,4,-1,21) ax.set_aspect (‘equal’) ax.set_title("Shear");Shear 4.4.5. Composition Represent the composition of two 2D transformations. Find the image of vector (10,0) when it is rotated by =/2 radians then stretched hori- zontally 2 units import numpy as np import matplotlib.pyplot as pit V = np.array([[2,311) origin = np.array([[0, 0, 0), (0, A=np.matrix((({0,-1],{1,0]]) Benp.matrix({(2,0],[0,1]]) Vienp.matrix(V) V2=Aenp. transpose (V1) 011) # origin point V2=np. array (V2) V3=np. array (V3) print ("Image of given vectors is:", V3) plt.quiver(+origin, V{:,0], VC:,1], color=['b'], scale=20) plt.quiver(+origin, V2(0,-), V2it,:], color=['r'], scale=20) plt.quiver(*origin, V3(0,:), V3(t,:], color=['g'], scale=20) plt.title('Blue=original, Red=Rotated ,Green=Rotated+Streached') plt.show()Blue=original, Red=Rotated ,Green=Rotated+Streached 0.08 002 0.08 “00s 00202 Another example. ¢ = np. array([[-cos (theta), sin(theta)],(sin(theta) , cos (theta)!1) C_coords = C@coords LTS © C_coords [0,:] yLLTS = Cleoords(1,:] # Create the figure and axes objects fig, ax = plt.subplots() # Plot the points. x and y are original vectors, x_LT1 and y_LT1 are dmages ax.plot (x,y, 'ro') ax. plot (x_LTS ,y_LTS, bo") # Connect the points by lines ax.plot (x,y, 'r',1s="--") ax. plot (x_LTS,y_LT8,'b') # Edit some settings ax. axvline(x-0,coler="k",18=":") ax. axhline(y*0, color="k",18=":") ax. grid (True) ax.axis((-2,2,-1,2]) ax.set_aspect (' equal’)4.5 Exercise: 1. Verify the rank nullity theorem for the following linear transformation a) T:R? + RS defined by T(z, y) = (x + dy, 2x + 5y, 3x + 6y). Ans: Rank=2, Nullity=1, RNT verified b) TR + Rt defined by T(x, y, 2) = (n+4y—z, 20+5y +82, 3r+y+22, 2+y+2) Ans: Rank=3, Nullity=1, RNT verified 2. Find the dimension of the subspace spanned following set of vectors a) S=(1,2,3,4), (2,4,6,8), (1, 1,1,1) Ans: Dimension of subspace is 2 b) $= (1,~1,3,4), (2,1,6,8), (1,1,1,1), (3,3,3,3) Ans: Dimension of subspace is 3 3. Find the image of (1,3) under following 2D transformations a) Horizontal stretch b) Reflection ©) Shear d) RotationLAB 5: Computing the inner product and orthogo- nality 5.1 Objectives: Use python 1. to compute the inner product of two vectors. 2. to check whether the given vectors are orthogonal. 5.2 Inner Product of two vectors Find the inner product of the vectors (2, 1, 5,4) and (3, 4,7, 8) import mumpy as ap #initialize arrays A= np.array((2, 1, 8, 4]) B= np.array((3, 4, 7, 81) dot product output = np.dot (A, B) print (output) 7 5.3. Checking orthogonality Verify whether the following vectors (2, 1,5,4) and (3,4,7,8) are orthogonal, import numpy as np #initialize arrays A= ap-array([2, 1, 5, 4]) B= np.array([3, 4, 7, 8]) fdot product output = np.dot(A, 8) print('Inner product is :*, output) Lf output==0 print (‘given vectors are orthognal else print ('given vectors are not orthognal ') Inner product is : 77 given vectors are not orthognal5.4 Exercise: 1. Find the inner product of (1,2, 3) and (3, 4,5). Ans: 26 2. Find the inner product of (1, ~1,2,1) and (4,2, 1,0). Ans: 4 3. Check whether the following vectors are orthogonal or not a) (1,1,—1) and (2,3,5). Ans: True b) (1,0,2,0) and (4,2, ~2,5). Ans: True ©) (1,2,3,4) and (2,3,4,5) . Ans: FalseLAB 6: Solution of algebraic and transcendental equa- tion by Regula-Falsi and Newton-Raphson method 6.1 Objectives: Use python 1. to solve algebraic and transcendental equation by Regula-Falsi method. 2. to solve algebraic and transcendental equation by Newton-Raphson method. 6.2 Regula-Falsi method to solve a transcendental equation Obtain a root of the equation 2° — 2x — 5 = 0 between 2 and 3 by regula-falsi method. Perform 5 iterations. # Regula Falsi method from sympy import + ymbol('x') input('Enter the function ') #ix"3-2ex-5; function ambdify (x,g) asfloat(input('Enter a valus :')) #2 float(input('Enter b valus :')) #3 Gapus (Enter number of iterations :')) #5 for i in range(1,W+1) e= (art (b) “bet (a))/(#()-£(a)) 34 (CE (a) #4 (6) <0)) bee else print ('itration id the root #0.3f \t function value 0.3f \n'h Gye, £695 Enter the function x*#3-2"x-5 Enter a valus :2 Enter b valus Enter number of iterations :5 itration 1 ‘the root 2.059 function value -0.391 itration 2 the root 2.081 function value -0.147 itration 3 ‘the root 2.090 function value -0.055 itration 4 ‘the root 2.093 function value -2.020 itration 5 the root 2.094 function value -0.007 Using tolerance value we can write the same program as follows: Obtain a root of the equation x — 2x —5 = 0 between 2 and 3 by regula-falsi method. Correct to 3 decimal places# Regula Falei method while loop2 from sympy import * x=Symbol ('x') g =input('Enter the function ') #hx3-2¥x-5; function f=lambdify(x,g) float (input (‘Enter a valus :')) #2 befloat(input('Enter b valus :')) #3 N=fleat (input (‘Enter telarence :')) # 0.00t while (abs(x-¢)>=N) c= (ast (b) bet (a))/(£(b)-£(a))) 5 Af (CE (a) #f(€)<0)) bee else “4 print('itration iid \t the root %0.3f \t function value 70.3f \n Gye, £60); print('final value of the root is 10.5f'Mc) Enter the function xt*3-2*x-5 Enter a valus :2 Enter b valus :3 Enter tolarence :0.001 itration 1 the root 2.059 function value -0.391 itration 2 ‘the root 2.081 function value -0.147 itration 3 the root 2.090 function value -.055 itration 4 ‘the root 2.093 function value -0.@2¢ itration 5 ‘the root 2.094 function value -0.007 itration 6 ‘the root 2.094 function value -0.003 final value of the root is 2.09431 6.3. Newton-Raphson method to solve a transcendental equa- tion Find a root of the equation 32 = 5 iterations $2 +1, near 1, by Newton Raphson method. Perform from sympy import + x=Symbol ("x") g =input('Enter the function ') #/3x-cos(x)-1; function fo lambdify(x,g) jag = aife(g);af=Lambdify (x, dg) xO= float (input ("Enter the intial approximation ')); # x0=1 n= int(input('Enter the number of iterations ')); #n=8; for i in range(1,n+1) x1 = (x0 ~ (£ (x0) /at (x0))) print('itration %d@ \t the root %0.3f \t function value ¥0.3f \n Gi, xt,£G1))); #print all iteration value xO = xl Enter the function 3tx-cos(x)-1 Enter the intial approximation 1 Enter the number of iterations 5 itretion 1 the root @.620 function value 0.046 itration 2 the root @.607 function value 0.008 itration 3 the root @.607 function value @.080 itration 4 the root @.6@7 function value 0.000 itration 5 the root @.607 function value 0.008 6.4 Exercise: 1. Find a root of the equation 32 = cos. +1, between 0 and 1, by Regula-falsi method. Perform 5 iterations. Ans: 0.607 Find a root of the equation ce* = 2, between 0 and 1, by Regula-falsi method. Correct to 3 decimal places. Ans: 0.853 Obtain a real positive root of x6 — x Perform 4 iterations. Ans: 1.856 0, near 1, hy Newton-Raphson method. Obtain a real positive root of 2*-+.2°—7x? 245 = 0, near 3, by Newton-Raphson method. Perform 7 iterations Ans: 2.061LAB 7: Interpolation /Extrapolation using Newton’s forward and backward difference formula 7.1 Objectives: Use python 1. to interpolate using Newton’s Forward interpolation method. 2. to interpolate using Newton’s backward interpolation method. 3. to extrapolate using Newton's backward interpolation method 1. Use Newtons forward interpolation to obtain the interpolating polynomial and hence 13 °5 7 9 calculate y(2) for the following: 6 10 62 210 S02 from sympy import + import numpy as ap n= int(input('Enter number of data points: ')) 210 x = np. zeros ((n)) y = np.zeros((n,n)) # Reading data points print (‘Enter data for x and y: ‘) for i in range(n) x[i] = float (input ( txCt+strG)e']=")) ylil [0] = floatCinpur( ty['+str(i)+"]=9)) # Generating forvard difference table for i in range(1,n) for j in range(0,n-i) y(gl(i) = yljrt) G1) ~ yi G10 print ('\nFORWARD DIFFERENCE TABLE\n'); for i in range(0,n) print ('0.2f" %(x[il), em for j in range(0, n-i) print ('\t\tno.2f" %¢y(a)(j)), end= print Q # obtaining the polynomial t=symbols('t') f-(] #f is a list type data ny p=(t-x(0})/(x(1)-x10]) F append (p) for {in range (1,o-1) f append (# [121] #(p-19/ (419) poly-y(0) [0] for {in range(o-1) | poay=poryey ol [aon] af (4 80sinp_poly-simplify (poly) print ('\aTHE INTERPOLATING POLYNOMIAL IS\n") ; pprint (simp_poly) # if you want to interpolate at sone point the next session will help imter=inpus("Do you want to interpolate at a point(y/a)? ") # y Mf inter=='y! a=float (input (‘enter the point ')) #2 interpol=lambdify (t, simp_poly) result» interpol (a) print('\nThe value of the function at! ,a,'is\n',result); Enter number of data points: § Enter data for x and yt FORWARD DIFFERENCE TABLE 1.00 6.00 4.00 48.00 48.00 0.00 3.00 10.00 32.00 96.00 28.00 5.08 62.00 148.08 144,00 7.00 210.00 292.00 9.00 502.00 THE INTERPOLATING POLYNOMIAL 1S 3 2 Lot -3.0t +10 + 7.0 Do you want to interpolate at a point(y/n)? y enter the point 2 The value of the function at 2.0 is 5. 2. Use Newtons backward interpolation to obtain the interpolating polynomial and x 135 7 9 hence calculate y(8) for the following data: G1 2 20 502 from sympy import * import mumpy as ap import sys print("This will use Newton's backvord intepolation formula ") # Reading number of unknowns n= intCinput('Enter number of data points: ')) # Making numpy array of n & n x n size and initializing # to zero for storing x and y value along with differences of y x = np. zeros ((n)) y = mp. zeros ((n,n)) # Reading data points 81print('Enter data for x and y: ') for i in range(n) x{i] = fleat(input( *x('+str(i)+']=")) yli} [0] = float(input( 'y('=str(a)+']=")) # Generating backward difference table for i in range(1,n) for j in range(m-1,i-2,-1) yCjlCil © y(g)la-a) = yCp-a) Ca-a] print ('\nBACKWARD DIFFERENCE TABLE\n'); for 4 in range(O,a) print (0.28 %G(4]), end for j in range(0, iv1) print \eno.26" A¢y Ci print, # obtaining the polynomial t=symbois('t') fe pe (tox (n-11)/ Celt] x10] ) f-append(p) for i in range(1,n-1) f. append (# [4-1] *(p+a)/(a+1)) poly-y(n-1] [0] print (poly) for i in range(n-1) poly=poly+y(n-1] [isi] +f (41 simp_poly=sinplify(pely) print \nTHE INTERPOLATING POLYNOMIAL 1S\n'); pprint (simp_poly) # if you want to interpolate at some point the next session will help inter=inpus('De you vant to interpolate at a point(y/n)? ') if inter=='y' afloat (input (‘enter the point ')) interpol=lambéity (t, simp_poly) result=interpol (a) print('\nThe value of the function at’ ,a,'is\n',result); a2This will use newton": backward intepolation formula Enter number oF data points: 5 inter data for x andy! x{e)-2 yie}=s x11)3 xp)-s yl2}=s2 ypj-210 xta)=9 yla}=502 BACKWARD DIFFERENCE TABLE 6.09 10.00 4.00 62.00 52.00 48.00 210,08 148,00 96.00 48,09 592,00 292.00 144,00 48.00 0.00 THE INTERPOLATING POLYHOMIAL 15 3 2 Let -3.0t +1004 7.0 0 you want to interpolate at 9 point(y/n)? y enter the point & The value of the function at 8.0 Is 335.0 7.2 Exercise: 1. Obtain the interpolating polynomial for the following data x 0 123 yo 121 10 Ans: 2x9 — 72? +60 +1 2. Find the number of men getting wage Rs. 100 from the following table wage: 50 150 250 350 No. ofmen: 9 30 35 42 Ans: 23 men 3. Using Newton’s backward interpolation method obtain y(160) for the following data x: 100 150 200 250 300 y: 1 13 15 17 1B Ans: 13.42 4, Using Newtons forward interpolation polynomial and calculate y(1) and y(10) x: 3 4 5 G6 7 8 9 y: 48 84 145 236 362 528 73.9 Ans: 3.1 and 100 83LAB 8: Computation of area under the curve using Trapezoidal, Simpson’s Qy" and Simpsons (3) rule 8.1 Objectives: Use python 1. to find area under the curve represented by a given function using Trapezoidal rule 5 . a ay 2. to find area under the curve represented by a given function using Simpson’s (4) 3 rule. 3. to find area under the curve represented by a given function using Simpson’s (3)"" rule. 4. to find the area below the curve when discrete points on the curve are given. 8.2 Trapezoidal Rule 5 Evaluate [phy a # Definition of the function to integrate ef my_func (x) return 1 / (1 + x #* 2) # Function to implement trapezoidal method def trapezoidal (x0, xn, n) he Gm = x0) fa # Calculating step size # Finding sum integration = my_func(x0) * my_func(xn) # Adding first and last terns for i in range(1, n) ke xO+ith # i-th step value Antegration = integration + 2 my_func(k) # Adding areas of the trapezoids # Proportioning sum of trapezoid areas integration = integration + h / 2 return integration # Input section lower_limit = float (input("Enter lower limit of integration: ")) upper_linit = float (input ("Enter upper limit of integration: ")) sub_interval = int(input("Enter number of sub intervals: ")) # Call trapezoidal() method and get result result = trapezoidal(lower limit, upper_limit, sub_interval) # Print result print ("Integration result by Trapezoidal method is: " , result) 84Enter lover Limit of integration: @ inter upper Linit of integration: 5 Enter number oF sub intervals: 10 Integration result by Trapezoidal sethod is: 1.3731¢80612301099 8.3. Simpson’s (3)"" Rule 5 Evaluate [ pbs # Definition of the function to integrate ef my_func (x) return 1 / G1 + x #* 2) # Function to implement the Simpson's one-third rule def simpson13(x0,xn,n) n= (xn = x0) fa # calculating step size # Finding sum integration (my_fune (x0) + my_func(xn)) k= x0 for i in range(1,a) if a2 == 0 integration = integration + 4 + my_func(k) else integration = integration + 2 * my_tunc(k) Roh # Finding final integration value integration = integration = h = (1/3) return integration # Input section Lower limit = float (input ("Enter lower limit of integration: ")) upper_linit = float (input ("Enter upper limit of integration: ")) sub_interval = int(input("Enter number of sub intervals: ")) # Call trapezoidal () method and get result result = simpson13(lower_limit, upper_limit, sub_interval) print("Integration result by Simpson's 1/3 method is: /0.6f" % (result) ) enter lover Limit of integration: 0 enter upper Limit of integration: 5 Enter nunber of sub intervals: 100 Integration result by Simpson's 1/3 method is: 1.408120 8.4 Simpson’s 3/8th rule Evaluate fj) >Lrdr using Simpson's 3/8 th rule, taking 6 sub intervals def simpsons_3_8_rule(f, a, b, n): 85(b- ad /n f(a) + £00) for i in range(1, m, 3) t t s+ 3 tat ish) ri in range(3, a-1, 3) ste 3" fara b) ri in range(2, a-2, 3) st52* tata hd returns * 342 / 8 get (x) return 1/(1tx*#2) # function here # lower limit upper limit ° 6H 6 # number of sub intervals result = simpsons.3_8_rule(f, a, b, n) print (3. 5f'/iresult) 1.27631 8.5 Exercise: 1 1. Evaluate the integral [ as using Simpson’s 2 rue luate the integral f der using Simpson's 5 a Ans: 0.23108 3 oe 2. Use Simpson's 5 rule to find / ix by taking seven ordinates, a Ans: 0.5351 3. Evaluate using trapezoidal rule | sintadr. Take n = 6 Ans: 1/2 4. A solid of revolution is formed by rotating about the 2-axis, the area between the x-axis, the lines x = 0 and x = 1, and a curve through the points with the following co-ordinates: x y 0.00 | 1.0000 0.25 | 0.9896 0.50 | 0.9589 0.75 | 0.9089 1.00 | 0.84151 Estimate the volume of the solid formed using Simpson’s 3rd rule. Hint: Required volume is fy? * dx. **[Ans: 2.8192|** 5. The velocity v(km/min) of a moped which starts from rest, is given at fixed intervals of time t(min) as follows. t 2 4 6 8 10 12 M 16 18 20 v: 10 18 25 29 32 2 1 5 2 0 Estimate approximately the distance covered in twenty minutes: Answer for 5. We know that ds/dt=v. So to get distance (s) we have to integrate. Here h = 2.2, v9 =0, v1 = 10, vp = 18, v3 = 25 ete # ve shall use simpson‘s 1/3 rule directly to estimate nea (0, 10 ,18, 28, 29,32 ,20, 11 ,8 2, 01 result=(h/3)*((y [0] +y (10) )+4=(y(1]+y (3]+y (5) +y [7] +y[9))+2*(y (2) ~y (4) +y [ el+y(81)) print (3. $f! result, "km.") 309.33333 km a7LAB 9: Solution of ODE of first order and first degree by Taylor’s series and Modified Euler’s method 9.1 Objectives: Use python 1. to solve ODE by Taylor series method. 2. to solve ODE by Modified Euler method, 3. to trace the solution curves, 9.2 Taylor series method to solve ODE Solve: % — 2y — 3e® with y(0) — 0 using Taylor series method at x — 0.1(0.1)0.3. ## module taylor 'OX,Y = taylor (deriv, x,y, xStop,h) 4th-order Taylor series method for solving the initial value problem {y 3} = (F(x, {y})}, where fy} = fylol,yl1],...yfn-1]> x,y = initial conditions xStop = terminal value of x = increment of x from aumpy import array aof taylor (deriv,x,y,xStop,n) x= (] y= X. append (x) Y.append(y) ubile x < xStop # Loop over integration steps D = deriv(x,y) # Derivatives of y We 1.0 for j in range (3) ¥ Build Taylor series A= Hen/(j + 1) ys y * DUjl*H #H = hoj/j! xox th X.append(x) # Append results to Y.append(y) # lists X and ¥ return array(X),array(¥) # Convert lists into arrays # deriv = user-supplied function that returns derivatives in the 4x fy'tol y'C1] y'f27 yi (n-1] yiCo) y'' C4) y' 02]... y'\fn-2) yi (ol yi C2) y! (2) yi''[n-1] yeeecoy yt yi 2) .. yt tnt) def deriv(x,y) D = zeros((4,1)) 88[0] = [2*y[0] + 3*exp(x)) DI1] = [4*yl0]+ 9xexp(x)) D[2] = [B4y[0]* 21+exp (x) D(3] = [16¥y(0]+ 48xexp(x)) return D x= 0.0 # Initial value of x xStop = 0.3 # last value y = array([0.0]) # Initial values of y on # Step size » X,Y = taylor (deriv,x,y,xStop,h) print("The required values are :at x= //0.2f, ye/0.8f, x=/0.2f, y=/0.82, UO.28, y=H0.5f,x 2f, yeho (x60) ,¥(0) ,x(1) ,¥(1] , x12) ,¥(2) +03] ,¥(3) )) The required values are :at x= 0.00, x = 0.20, y=0.81079,x = 0.30, y=1.41590 0.00000, x=0.10, y=0.34850, Solve y’ | dy = 2? with initial conditions y(0) = 1 using Taylor series method at x = 0.1,0.2 from numpy import array def taylor (deriv ,x,y,xStop ,h) x= 01 yet) X. append (x) Y. append (y) while x < xStop # Loop over integration steps D = deriv(x,y) # Derivatives of y He 1.0 for j in range(3) # Build Taylor series R= Heas(j + 1) yoy + DUjlH # oH = h-g/j! eek th X.append(x) # Append results to Y.append(y) # lists X and Y return array (X),array(¥) # Convert lists into arrays # deriv = user-supplied function that returns derivatives in the 4 xn array fy‘ Co} y'l1) y' C2) yitn-t] y"Co] y"(1] y"(2) y"(n-1] yr'tol yer] yt 02]... yt fant] vento) yeeta) y*" C2) yeetn-t}) def deriv(x,y) D = zeros ((4,1)) D(0) = (x*+2-4+y(0]) DUA] = [2ex-4ex++2+16+y C0 D[2] = [2-Bexti6exee2-64ey (01) D(B] = [-8+32%x-84x0=2+256%y (0)] 89return D x= 0.0 # Initial value of x xStop = 0.2 # last value y = array((1.01) # Initial values of y ho. # Step size X,Y = taylor (deriv,x,y,xStop,h) print ("The required values are :at x= /0.2f, y=/0.5£, x=/0.2£, y=/0.5£, x = (0.2£, y=K0.5£"/(XC0] ,¥ [0] x4 D.YE1) ,x(2),¥(21)) The required values are :at x= 0.00, y=1.00000, x=0.10, y=0.66967, x = 0.20, y=0.45026 9.3. Euler’s method to solve OD ‘To solve the ODE of the form # = f(x,y) with initial conditions y(xo) = yo. The iterative formula is given by : y(eqsay = ylai) + A(x, y(2s)) Solve: y’ = e~* with y(0) = —1 using Buler’s method at x = 0.2(0.2)0.6. import numpy as mp import matplotlib.pyplot as plt # Define parameters f= lambda x, y: mp.exp(-x) # ODE h = 0.2 # Step size yO = -1 # Initial Condition ao3 # Explicit Euler Method yo] = yo x{0]+0 for i in range(0, n) xlitt]=x(il+h yli +a] = y[a] + met(xtal, ylal) print("The required values are at x= /0.2f, y x = 10.2f, yeK0.8f,x = 10.2, yaio 5£%(x(0] ,y [0] ,x{1],y[1],x{2] ,y(2], x(3),y(3])) prin("\n\n") plt.plot(x, y, 'bo--', label='Approxinate') plt.ploc(x, -np.exp(-x), ‘g*-', label='Exact') plt.titie ("Approximate and Exact Solution" ) plt.xlabel('x') plt.ylabel('£(x)') plt.grid() plt.logend(1oc=' best") plt.shov()The required values are at x= 0.00, y= x = 0.40, y=-0.63625,x = 0.60, y=-0.50219 20, y=-0.80000, Approximate and Exact Solution ~e- Approximate 35 fe Bxact 25 20 tbo 1 10 9) 0 20 20 60 80 100 Solve: y’ = —2y + 2%e~™ with y(0) = 1 using Buler’s method at x = 0.1,0.2 import numpy as np import matplotlib.pyplot as plt # Define parameters £ = lambda x, y: ~2+y+(x+*3)=np.exp(-2*x) # ODE h = 0.1 # Step size yO = 1 # Initial Condition * 2 Explicit Euler Method y(0] = yo x(0]=0 for i in range(0, n) xlattl=x(i]+h yli + 4] = yl] + nef(xCa, yCal) print("The required values are at x= /0.2f, 0. 5£\n\n"%(x[0] ,y [0] x12 dsy(t),x(21,y (21) plt.plot(x, y, ‘bo--', label="Approximate (Euler's method)") plt.title("Solution by Euler's method") plt.xlabel('x') plt.ylabel('£(x)") plt-gridQ) plt. legend(1oc='best') plt. show ()The required values are at x= 0.00, 20.20, y=0.64008 00000, x=0.10, y=0.80000, Solution by Euler's method e- Approximate(Euler's method) a 3 30 25 20 fix) 3 10 0 20 40 60 20 100 9.4 Modified Euler’s method ‘The iterative formula is. os h ; i = yo FU (em) + Flea vi Ih m= 041, 2.3.--04 where, y{” is the n'* approximation to yy ‘The first iteration will use Euler's method: y!°! = yo + hf (0, yo). Solve y’ = —ky with y(0) = 100 using modified Euler’s method at 2 = 2 100, by taking import numpy as np import matplotlib.pyplot as plt def modified euler(f, x0, yO, hy n) x = np.zeros(n*1) y = mp.zeros (art) x [0] = x0 ylo] = yo i in range(n) xlitt] = xi] +h kl © bh * f(x(i], y(il) k2 = ho * £(x[i+i], yli) + kt) ylitt] = yli) + 0.8 = Ct + 42) return x, ydet f(x, y) return -0.01 + y # ODE dy/ax = ~Ky x0 = 0.0 yo = 100.0 25 4 print("The required value at x= 10.28, y=/0.5£"%(x[4) ,y(4])) print("\n\n") # Plotting the results pit.plov(x, y, 'bo-') pit .xlabel ("x") pit. ylabel('y') pit.title('Solution of dy/dx = -ky using Modified Euler\'s Method") plt.grid(True) |plt-shovO bh x, y = modified euler(f, x0, yO, h, a) ‘The required value at x= 100.00, y=37.25290 Solution of dy/dx ky using Modified Euler's Method 100 80 70 60 o 20 0 60 80 300 9.5 Exercise: 1. Find y(0.1) by Taylor Series exapnsion when y' Ans: y(0.1) = 0.9138 xv y(0) ) ) 2. Find y(0.2) by Taylor Series exapnsion when y! = 2?y — 1,y(0) =1,h =0.1 Ans: (0.2) = 0.802273. Evaluate by modified Buler’s method: y/ = In(x + y), y(0) = 2 at x = 0(0.2)0.8. Ans: 2.0656, 2.1416, 2.2272, 2.3217 4. Solve by modified Euler’s method: y/ = 2+ y,y(0) = 1,h = 0.1,2 = 0(0.1)0.3. Ans: 1.1105, 1.2432, 1.4004 4LAB 10: Solution of ODE of first order and first de- gree by Runge-Kutta 4th order method and Milne’s predictor and corrector method 10.1 Objectives: 1. To write a python program to solve first order differential equation using 4th order Runge Kutta method 2. To write a python program to solve first order differential equation using Milne’s predictor and corrector method. 10.2 Runge-Kutta method Apply the Runge Kutta method to find the solution of dy/dx = 1+ (y/2) at y(2) taking h— 0.2. Given that y(1) = 2. from sympy import * import numpy as np def RungeKutta(g,x0,h,y0, xn): x,yssymbols('x,y") f-lambdity((x,y],¢) xt=x0+h Y= [yl while xt<=xn wish4f (x0, 0) R2-ReE(KO+R/2, yOKL/2) K3=h+E(xO+h/2, yO+k2/2) wAshst(xOrh, yOPKS) yl-yO1 (1/6) + (kL 2+k2 +243 +4) Y. append (y1) aprint ('y (3. 3f"ixt,") is U3. 3f"hyt) x0=xt yoryt xextth return np.round(¥,2) RungeKutta('1+(y/x)',1,0.2,2,2) array((2., 2.62, 3.27, 3.95, 4.66, 5.39]) 10.3. Milne’s predictor and corrector method Apply Milne’s predictor and corrector method to solve dy/dr = x? + (y/2) at y(1-4) Given that y(1)=2, y(1.1)=2.2156, y(1.2)=2.4649, y(1.3)=2.7514. Use corrector formula, thrice. # Milne's mothod to solve first order DE # Use corrector formula thrice y' 95yi-2.2156 ya=2. 4649 y3=2.7514 h-o.1 x1=x0¢h xQexi +h x3=x2¢h xé=x3+h def f(x,y) return xe+2+(y/2) y0-£(x0, 0) yitet Ga, yt) (2,92) (3,93) yap~y0" (4h/3) +(2¥y11-y1202«y13) print('predicted value of y4 is 13.3£"/y4p) y14=£ G4, yA); for i in range (1,4) y4oy2+(n/3) + (ylaranyt3ey42) 5 print (‘corrected value of y4 after \t iteration Zid is \t 13.5f\t Gay) yla=t (x4, y4); predicted value of y4 is 3.079 corrected value of y4 after iteration 1 is 3.07940 corrected value of y4 after iteration 2 is 3.07940 corrected value of y4 after iteration 3 is 3.07940 In the next program, function will take all the inputs from the user and display the answer Apply Milne’s predictor and corrector method to solve dy/de = x? + (y/2) at y(1.4) Given that y(1)=2, y(1.1)=2.2156, y(1.2)=2.4649, y(1.3)=2.7514. Use corrector formula, thrice. ‘von sympy import + def Milne(g,x0,h.y0,y1.y2+¥3) x,y-symbois (x, 9") flambdify ([x,y] 46) xi=x0+h y10=£(x0, yo) yil-£G1,y1) yl2=£ (32, y2) y13=£ (3, y3) yap~yO! (44h / 3) (Qeyit-yt2+2+y13) print('predicted value of y4",y4p) yl4=t (x4, yap) for i in range (1,4) yAey2+(h/3) (y14eaey13-y12) print (corrected value of y4 , iteration a ‘Ki,y4 96yla=t (x4, y4) Milne (‘x4*2y/2',1,0.1,2,2. 2156 ,2. 4649, 2.7514) predicted value of y4 3.0792733333333335 corrected value of y4, iteration 1 3.0793962222222224 corrected value of y4, iteration 2 3.079398270370371 corrected value of y4 , iteration 3 3.079398304506173 Apply Milne’s predictor and corrector method to solve dy/dx = x — y* , y(0)=2 obtain y(0.8). Take h=0.2, Use Runge-Kutta method to calculate required initial values YoRungeKutta('x-y+#2',0,0.2,0,0.8) print('y values from Runge ~Kutta method: ',Y) Mitne('x-y##2',0,0.2,¥(0) ,¥(1),¥(2),¥(3]) y values from Runge -Kutta method: [0. 0.02 0.08 0.18 0.3] predicted value of y4 0.3042133333333334 corrected value of y4 , iteration 1 0.3047636165214815 corrected value of y4 , iteration 2 0.3047412758696499 corrected value of y4 , iteration 3 0.3047421836520892 10.4 Exercise: 1. Find y(0.1) by Runge Kutta method when y/ = x — y?, y(0) Ans: y(0.1) = 0.91379 2. Evaluate by Runge Kutta method : y/ =log(x + y),y(0) = 2 at x = 0(0.2)0.8 Ans: 2.155, 2.3418, 2.557, 2.801 3. Solve by Milnes method: y’ = x+y, y(0)=1, h=0.1, Calculate y(0.4) . Calculate required initial values from Runge Kutta method. Ans: 1.583649219