Python for Rapid Engineering Solutions - Finding Roots - PDF

Document Details

DarlingBowenite8037

Uploaded by DarlingBowenite8037

Arizona State University

Steve Millman

Tags

python programming numerical methods root finding optimization

Summary

This document is lecture notes on using python for rapid engineering solutions. It focuses on the topic of finding roots using numerical methods, specifically Newton's Method and fsolve. It provides code examples and explanations to help learners quickly grasp the topic. It further introduces the Relaxation method and concepts for optimizations.

Full Transcript

Python for Rapid Engineering Solutions Steve Millman Finding Roots Welcome! Today’s Objectives: Finding roots of nonlinear equations Examine numerical methods Use fsolve to find roots Look at some examples Example code from Newman Solving Nonlinear Equations Many equations of...

Python for Rapid Engineering Solutions Steve Millman Finding Roots Welcome! Today’s Objectives: Finding roots of nonlinear equations Examine numerical methods Use fsolve to find roots Look at some examples Example code from Newman Solving Nonlinear Equations Many equations of interest require numerical methods Newton's Method: 𝑓(𝑥) 𝑓! 𝑥 = ∆𝑥 𝑓(𝑥) ∆𝑥 = 𝑓′(𝑥) 𝑓(𝑥) 𝑓(𝑥) 𝑥 − 𝑥"#$ = 𝑓′(𝑥) 𝑓(𝑥) 𝑥𝑛𝑒𝑤 𝑥 𝑥𝑛𝑒𝑤 = 𝑥 − root 𝑓′(𝑥) Newton's Method Advantage – it converges quickly Disadvantages: Have to know 𝑓 ) 𝑥 There may be multiple roots It jumps too far when 𝑓 ) 𝑥 is small Possible Solution: add a factor r to limit ∆𝑥 𝑓(𝑥) 𝑥𝑛𝑒𝑤 = 𝑥 − 𝑟 𝑓′(𝑥) Newton's Method Code: Set Up m8_newton.py: import numpy as np # import math and arrays import matplotlib.pyplot as plt # so we can plot from scipy.optimize import newton # newton method function ################################################################################ # Function which has a zero when x is an odd multiple of pi # # input: x; returns computed value # ################################################################################ def funcpi(x): return(1.0+np.cos(x)) ################################################################################ # Function which is the derivative of funcpi # # input: x; returns the derivative of funcpi at x # ################################################################################ def dfuncpi(x): return(-np.sin(x)) Newton's Method Code: Local Version m8_newton.py (continued): ################################################################################ # A user-defined Newton's Method # # Inputs: # # func: the function to use # # dfunc: its derivative # # xstart: starting value, that is, the first guess # # rconv: convergence rate # # Outputs: # # xn: the root # # num: number of iterations # # # # NOTE: There is no safety in this function in case of nonconvergence! # ################################################################################ def newtmethod(func,dfunc,xstart,rconv,tol): xn = xstart # initialize the value delta = 1 # initialize the difference num = 0 # initialize the iterations while(delta > tol): # while too big... xnp1 = xn - rconv*func(xn)/dfunc(xn) # compute next value delta = abs(xnp1 - xn) # compute the differnce xn = xnp1 # swap to the next guess num += 1 # increment the counter return(xn,num) # return the results Newton's Method Code: Call the Function m8_newton.py (continued): # Find the zero for the function # sensitive to r! start at 0.1, convergence parameter 1 soln,num = newtmethod(funcpi,dfuncpi,0.1,1,1e-6) print(' estimate pi ', soln, num) # this time, start at 0.1, but convergence parameter 0.1 soln,num = newtmethod(funcpi,dfuncpi,0.1,0.1,1e-6) print(' estimate pi ', soln, num) # Now use the built in newton's method soln_ = newton(funcpi,0.1,dfuncpi,tol=1e-6) print(" built in newton's method ", soln_) # plot the function X = np.arange(0,25,0.01) Y = funcpi(X) plt.plot(X,Y) plt.grid() plt.title("1+cos(x)") plt.show() Newton's Method Code: Run It Ø python m8_newton.py estimate pi 21.991147651841388 21 estimate pi 3.141573848820656 213 built in newton's method 21.991147651841388 Relaxation Method Rearrange the equation to get the form 𝑥 = 𝑓(𝑥) 1. Start with a guess for 𝑥 2. Plug the guess into 𝑓(𝑥) to get a new value of 𝑥 3. Repeat step 1 & 2 with the new guess until ∆𝑥 < 𝜀 Relaxation Method Example Example: 𝑥 + + 2𝑥 + 1 = 0 (roots are -1!) , ! -. Rewrite as 𝑥 = − + Disadvantages: Multiple roots Doesn't converge Can't get to 𝑥 = form Relaxation Code: Set Up m8_relax.py: START_VALS = [2,.5] # pick some starting values TOO_BIG = -1e35 # if answer gets bigger than this, stop! ITERATIONS = 20000 # maximum number of iterations REPORT_ALL = 5 # report all values for first few iterations REPORT_MOD = 3000 # then report only occasionally ################################################################################ # Compute a new value of x for the relaxation method # # input: x_in, the current value of x # # output: the new value of x # ################################################################################ def relax(x_in): return -((x_in * x_in) + 1)/2 Relaxation Code: Calculation m8_relax.py (continued): # keep iterating from the start value until done for x in START_VALS: print("Starting value",x) error = 0 for iter in range(ITERATIONS): x = relax(x) if ( iter < REPORT_ALL ) or ( iter % REPORT_MOD == 0 ): print(iter,":",x) if x < TOO_BIG: print(iter,":",x,"Not Converging!") error = 1 break; if not error: print("Final value after",ITERATIONS,"iterations",x,"\n") else: print("") Relaxation Code: Run It Ø python m8_relax.py Starting value 2 0 : -2.5 1 : -3.625 2 : -7.0703125 3 : -25.494659423828125 4 : -325.4888295684941 8 : -4.843510864363938e+35 Not Converging! Starting value 0.5 0 : -0.625 1 : -0.6953125 2 : -0.741729736328125 3 : -0.7750815008766949 4 : -0.800375666500635 3000 : -0.99933594045626 6000 : -0.9996673579800254 9000 : -0.9997780952088215 12000 : -0.9998335159345848 15000 : -0.9998667855353348 18000 : -0.9998889725717806 Final value after 20000 iterations -0.9999000633203341 Binary Search Find an 𝑥1 such that 𝑓 𝑥1 < 0 Find an 𝑥2 such that 𝑓 𝑥2 > 0 ,$-,% 1. Set 𝑥3 = + 2. if 𝑓 𝑥3 < 0 set 𝑥1 = 𝑥3 else set 𝑥2 = 𝑥3 𝑥1 3. Repeat until good enough 𝑥2 fsolve scipy.optimize.fsolve finds roots Don't confuse this with numpy.linalg.solve! fsolve(function,initial_guess,args=(arg_list)) function is the function whose roots are desired initial_guess may be an array arg_list are any additional arguments not being varied Example: Earth-Moon Lagrange Point Find the point between the earth and moon where their gravitational forces cancel. Moon Put satellite of mass msat at the Lagrange point R-r 𝐺𝑀𝑚%&' 𝐺𝑚𝑚%&' Earth mass m r ( − ( = 𝑚%&' 𝜔( 𝑟 𝑟 𝑅−𝑟 R 𝐺𝑀 𝐺𝑚 − − 𝜔( 𝑟 = 0 𝑟( 𝑅−𝑟 ( mass M Lagrange Point Code: Set Up m8_lagrange.py: from scipy.optimize import fsolve # root finder import matplotlib.pyplot as plt # plot results import numpy as np # for arrays G = 6.674e-11 # gravitational constant ME = 5.974e24 # mass of earth in kg MM = 7.348e22 # mass of moon in kg R = 3.844e8 # distance to moon in m OMEGA = 2.662e-6 # rad/s of moon ################################################################################ # Function which is force between earth and moon to find the zero # # Input: # # r - the radius to try, that is, distance from earth # # Output: # # func - the result of the function - trying to make it 0 # # # # Note that the equation used was derived in class... # ################################################################################ def lagrange_calc(r): func = ( ( G * ME ) / ( r**2 ) ) - \ ( ( G * MM ) / ( ( R - r )**2 ) ) - \ ( r * OMEGA**2 ) return func Lagrange Point Code: Solve and Graph m8_lagrange.py (continued): # now find the root using an initial guess of 3e8, which is ~78% to the moon root = fsolve(lagrange_calc,3e8) # Note that we know that ths problem has exactly one root... # (fsolve returns an array, but in this case it has only one entry.) print("distance from earth: {:,}".format(int(root)), "meters") print("distance from moon: {:,}".format(int(R-root)), "meters") # create a plot showing the function across most of the range x = np.arange(0.5e8, 3.8e8, 0.1e7) plt.plot(x, lagrange_calc(x)) plt.plot([root],,marker='x',markersize=20,c='red') plt.plot([R],,marker='o',markersize=5,c='black',label='distance to moon') plt.xlabel('distance from sun') plt.ylabel('residual') plt.xlabel('distance from earth') plt.ylabel('residual') plt.title('Earth-Moon Lagrange Point') plt.show() Ø python m8_lagrange.py distance from earth: 326,045,071 meters distance from moon: 58,354,928 meters Lagrange Point Code: Visualization Sun-Earth Lagrange Point A few simple changes of constants! Ø python m8_lagrange_sun.py distance from sun: 147,993,882,983 meters distance from earth: 1,506,117,016 meters Sun-Earth Calculation Lagrange Points This Photo by Unknown Author is licensed under CC BY-SA Python for Rapid Engineering Solutions Steve Millman Finding Roots and Minima and Maxima Welcome! Today’s Objectives: Another fsolve example Examine a method for finding minima and maxima Example code from Newman Diode Example: Solve a Circuit Find the voltages V1 and V2. VDD ' 𝑛𝑘𝑇 R1 R3 𝐼!"#!$ = 𝐼% 𝑒 &'! − 1 𝑤ℎ𝑒𝑟𝑒 𝑉( = ≈ 0.5 𝑞 V1 V2 Node currents sum to 0: R2 R4 𝑉) − 𝑉** 𝑉) + + 𝐼!"#!$ = 0 𝑅) 𝑅+ 𝑉+ − 𝑉** 𝑉+ + − 𝐼!"#!$ = 0 𝑅, 𝑅- Use fsolve to solve two simultaneous equations Diode Code m8_diode.py: from scipy.optimize import fsolve # root finder import numpy as np # for math VDD = 5 # supply voltage R1 = 1e3 # various resister values R2 = 4e3 Ø python m8_diode.py R3 = 3e3 [3.44695462 2.82956807] R4 = 2e3 I0 = 3e-9 # diode saturation current VT = 0.05 # diode threshold voltage ################################################################################ # fsolve will try to find values of v to drive these two equations to 0 # # Input: # # v - a list of two entries representing the voltages of the two nodes # # Output: # # a, b - the results of apply the voltages to the equations - want 0! # ################################################################################ def f(v): diode = I0*(np.exp((v-v)/VT)-1.0) # current through the diode a = (v-VDD)/R1 + v/R2 + diode # solve node 1 b = (v-VDD)/R3 + v/R4 - diode # solve node 2 return [a,b] print(fsolve(f,[4,2])) # returns the values for V1 and V2 in an array Finding Maximum: The Golden Search Select 4 points, T1 below, T4 above, and T2, T3 between If outer values bigger, then bad choices! Otherwise, squeeze! (Note: not the golden ratio search) If f(T2) > f(T3): T1 T2 T3 T4 Slide T3 and T4 down T1 T2 T3 T4 T2 splits the distance to T1 Else: T1 T2 T3 T4 Slide T1 and T2 up T1 T2 T3 T4 T3 splits the distance to T4 The Golden Search Code m8_golden.py: ################################################################################ # Function to implement the golden search # # Inputs: # # func - the function for which a max is desired # # x1,x2,x3,x4 - the current values, lowest to largest # # tol - the tolerance: quit if x1 and x4 are closer than tol # # Outputs: # # x1,x2,x3,x4 - the new values, lowest to largest # ################################################################################ def goldsearch(func,x1,x2,x3,x4,tol): if (x4-x1>tol): # not close yet? if (max(func(x2),func(x3))>max(func(x1),func(x4))): # a middle > outside if (func(x2)>func(x3)): # if x2 is bigger than x3... x4 = x3 # slide x4 down to x3 x3 = x2 # slide x3 down to x2 x2 = (x1+x3)/2.0 # x2 is average of x1 and old x2 else: # if x3 was bigger than x2... x1 = x2 # slide x1 up to x2 x2 = x3 # slide x2 up to x3 x3 = (x2+x4)/2.0 # x3 is average of x4 and old x3 # either way, we search again with the new values x1,x2,x3,x4 = goldsearch(func,x1,x2,x3,x4,tol) return(x1,x2,x3,x4) Calling the Golden Search m8_goldsimple.py: import numpy as np # for math and arrays from m8_golden import goldsearch # the golden search function ################################################################################ # A test function to find a max for... Clearly, the max is at x = 2. # # Input: # # x - value at which to evaluate the function # # Output: # # the value of the function evaluated at x # ################################################################################ def func(x): return(float(10.0-(x-2.0)**2)) x1,x2,x3,x4 = goldsearch(func,0,1,2,2.5,0.001) # now find the max print("max is at : ",max(x2,x3), "where it is f(x): ", func(max(x2,x3))) Ø python m8_goldsimple.py max is at : 2.0001220703125 where it is f(x): 9.999999985098839 Optimizing Light Bulb Temperature At what temperature is visible light maximized? 𝜆56 𝐼 𝜆 = 2𝜋𝐴ℎ𝑐 4 897 𝑒 :;! < −1 :# 𝜆56𝑑𝜆 , 897 :" 𝑒 :;! < −1 𝜂= @ 𝜆56𝑑𝜆 , 897 = >? 𝑒 :;! < −1 Simplify ℎ𝑐 ℎ𝑐 ℎ𝑐 𝐿𝑒𝑡 𝑥 = then λ = 𝑎𝑛𝑑 𝑑𝜆 = − 4 𝑑𝑥 𝜆𝑘A 𝑇 𝑥𝑘A 𝑇 𝑥 𝑘A 𝑇 6 B 56 𝑥𝑘A 𝑇 ℎ𝑐 𝑘A 𝑇 𝜆 𝑑𝜆 = × − 4 𝑑𝑥 = − 𝑥 C𝑑𝑥 ℎ𝑐 𝑥 𝑘A 𝑇 ℎ𝑐 Ignore the constants as they'll cancel, so denominator: @ @ 𝜆56𝑑𝜆 𝑥 C𝑑𝑥 𝜋B , 897 𝑏𝑒𝑐𝑜𝑚𝑒𝑠 , D = 𝑒 − 1 15 = >? 𝑒 :;! < −1 = >? Not So Lucky with the Numerator 897 :# :" ;! < 𝜆56𝑑𝜆 𝑥 C𝑑𝑥 , 897 𝑏𝑒𝑐𝑜𝑚𝑒𝑠 , 𝑒D − 1 :" 𝑒 :;! < −1 897 :# ;! < So need to solve: 897 :" ;! < 15 𝑥 C𝑑𝑥 , 𝜋B 𝑒D − 1 897 :# ;! < Bulb Problem: Set Up m8_bulb.py: import numpy as np # for math and arrays from scipy.integrate import quad # to integrate import matplotlib.pyplot as plt # to plot the curve from m8_golden import goldsearch # the golden search function # Constants used in the calculations LAMBDA1 = 430e-9 # lower wavelength LAMBDA2 = 750e-9 # upper wavelength BOLTZ = 1.38064852e-23 # Boltzmann Constant J/K C = 2.99792458e8 # Speed of light in vacuum m/s PLANCK = 6.626070040e-34 # Planck's Constant J s # these coefficients are used in the effincandescent_bulb function COEFF1 = PLANCK*C/(LAMBDA1*BOLTZ) # so only compute once! COEFF2 = PLANCK*C/(LAMBDA2*BOLTZ) # so only compute once! COEFF3 = 15.0 / ( np.pi**4.0 ) # so only compute once! ################################################################################ # Function to be integrated for lightbulb efficiency # # Input: x - variable inside integrand # # Output: function evaluated at x # ################################################################################ def expfun(x): return(x**3/(np.exp(x)-1.0)) Bulb Problem: Find the Maximum m8_bulb.py (continued): ################################################################################ # Function to perform the lightbulb efficiency integration # # Input: # # temp - temperature at which to do the calculation # # Output: # # the efficiency # ################################################################################ def effincandescent_bulb(temp): upperlimit = COEFF1/temp # calculate integration limits lowerlimit = COEFF2/temp res,err = quad(expfun,lowerlimit,upperlimit) # do the integration effic = COEFF3 * res # mult by constant out front return(effic) # find the temperature of bulb with the peak efficiency # first look at 300 K, where it is not very efficient temp = 300 eff = effincandescent_bulb(temp) print('efficiency ',eff, ' T ',temp) # now do the search T1,T2,T3,T4 = goldsearch(effincandescent_bulb,300,2500,7500,10000,0.001) Bulb Problem: Print Max and Actual m8_bulb.py (continued): T2_eff = effincandescent_bulb(T2) # efficiency at T2 T3_eff = effincandescent_bulb(T3) # efficiency at T3 if ( T2_eff > T3_eff ): peak_temp = T2 # this is the peak T and efficiency eff_peak = T2_eff * 100 else: peak_temp = T3 # this is the peak T and efficiency eff_peak = T3_eff * 100 print("max is at : "+ str(round(peak_temp,0)) + " where it is: " + str(round(eff_peak,1)) + '%') print('\nactual filament temperature 2000 to 3300 K') print('suggesting an efficiency from ' + str(round(100*effincandescent_bulb(2000),1)) + '% to '+ str(round(100*effincandescent_bulb(3300),1)) + '%' ) Bulb Problem: Visualize the Data m8_bulb.py (continued): # plot the efficiency over a range of temperatures from 300 to 10300K temps = np.linspace(300,10300,50,float) # create a range of interest Eff = [] # an empty list to hold efficiencies for a_temp in temps: # for each temperature... Eff.append(effincandescent_bulb(a_temp)) # append the efficiency at that T ax = plt.axes() ax.set_facecolor('grey') plt.title(' Efficiency vs Temperature ') # now create a plot plt.plot(temps,Eff) plt.ylabel(' Efficiency ') plt.xlabel(' Temperature K ') # place arrow at peak: plt.arrow(x,y,dx,dy) draws line plt.arrow(peak_temp,0,0,eff_peak,color='blue') # place arrow over range of operating temp of 2000 to 3300 eff_op_range = effincandescent_bulb(3300) plt.arrow(2000,0,1300,0,color='orange') plt.arrow(3300,0,0,eff_op_range,color='orange') plt.text(peak_temp+100,.2,"max efficiency",fontsize=24,color="blue") plt.text(3400,.1,"typical\nefficiency",fontsize=24,color="orange") plt.show() Bulb Problem: Run It Ø python m8_bulb.py efficiency 7.148e-22% at T = 300 max is at : 6561.0 where it is: 39.1% actual filament temperature 2000 to 3300 K suggesting an efficiency from 1.3% to 14.7% Python for Rapid Engineering Solutions Steve Millman Optimization Welcome! Today’s Objectives: Understand what optimization is Use different methods to optimize functions Optimization Goal: find the minimum or maximum of a function Methods from the last lecture were an introduction There are many, far more powerful, optimizers in Python scipy.org has an excellent chapter link to this on the course web site Problem: Optimization is very difficult! Find the minimum Easy to get trapped in local minima Many algorithms are available… The Six-Hump Problem: Set Up m8_sixhump.py: import numpy as np # for arrays from scipy import optimize # access to optimization functions import matplotlib.pyplot as plt # to create the plot from mpl_toolkits.mplot3d import Axes3D # to allow 3D! ################################################################################ # Function to implement the 6-hump shape from which to find the minimum # # Input: # # coords - x and y value at which to evaluate the function # # Output: # # the value of the function # ################################################################################ def sixhump(coords): x = coords # extract the x and y for clarity y = coords return ( 4 - ( 2.1 * x**2 ) + ( x**4 / 3. ) ) * x**2 + \ x * y + \ (-4 + ( 4 * y**2 ) ) * y**2 The Six-Hump Problem: 2-D Visualization m8_sixhump.py (continued): # create a meshgrid so we have two arrays from which we can project the surface x = np.linspace(-2, 2) # allow x to vary from -2 to 2 y = np.linspace(-1, 1) # allow y to vary from -1 to 1 xg, yg = np.meshgrid(x, y) # create the meshgrid plt.figure() # simple visualization in 2D plt.imshow(sixhump([xg, yg])) plt.colorbar() plt.xlabel('x') plt.ylabel('y',rotation=0) plt.title('six hump in 2D') plt.show() Ø python m8_sixhum.py The Six-Hump Problem: Optimize and 3-D m8_sixhump.py (continued): # First optimization using bfgs with initial guess 2,-1, disp=0 for quiet result = optimize.fmin_bfgs(sixhump, [2,-1], disp=0) print('bfgs result:',result) # Second optimization using basin hopping with initial guess 2,-1 result = optimize.basinhopping(sixhump, [2,-1]) print(type(result)) # this is a significant structure print('basin hopping result:',result.x) # many other attributes generated fig = plt.figure() # awesome visulization in 3D ax = fig.add_subplot(111, projection='3d') surf = ax.plot_surface(xg, yg, sixhump([xg, yg]), rstride=1, cstride=1, cmap=plt.cm.jet, linewidth=0, antialiased=False) ax.set_xlabel('x') ax.set_ylabel('y') ax.zaxis.set_rotate_label(False) ax.set_zlabel('f(x, y)') ax.set_title('Six-hump Camelback function') plt.show() bfgs result: [ 0.089842 -0.7126564] basin hopping result: [ 0.08984201 -0.71265641] The Six-Hump Problem: Graph Residuals: Minimizing the Difference Goal: Fit a curve to data; achieve a nonlinear fit! Solution: propose a function f and fine tune parameters Method: Use a least squares fit Required: A function that returns 0 if there's a match Implementation: leastsq calls the residual function the residual function calls the target function leastsq Example: Set Up Given data, can we determine kd? m8_leastsq.py: import numpy as np # math and arrays from scipy.optimize import leastsq # optimization function import matplotlib.pyplot as plt # so we can plot N=1000 # number of samples to generate ################################################################################ # Function to fit # # Inputs: # # kd - first value, which we are trying to determine # # p0 - second value, which is the known domain # # Output: # # value of the function # ################################################################################ def func(kd,p0): return 0.5*(-1-((2*p0)/kd) + np.sqrt(4*(p0/kd)+(((2*p0)/kd)-1)**2)) leastsq Example: The Residual Function m8_leastsq.py (continued): ################################################################################ # Function to compute the difference between the actual and predicted values # # Inputs: # # kd_guess - the guess for the value of the first parameter # # p0 - the second parameter, which is known # # actual - the sampled value # # Output: # # difference between the actual value and the value calculated with guess # ################################################################################ def residuals(kd_guess,p0,actual): return actual - func(kd_guess,p0) leastsq Example: Analysis m8_leastsq.py (continued): # Create a noisy signal based on a known value of kd of 3 # since random returns uniformly distributed [0,1), subtract.5 gives # [-.5,.5) so we get +/- noise kd=3.0 # the "known" value p0 = np.linspace(0,10,N) # an array for p0 clean = func(kd,p0) # the clean signal actual = clean+(np.random.random(N)-0.5)*1.0 # the noisy signal # now try to extract the known value of kd by minimizing the residuals # residuals - the function we are optimizing # 5 - the initial guess for the value to find, in this case, kd # args - the additional arguments needed, in this case the x and y values # full_output - return all the outputs kd_match,cov,infodict,mesg,ier = leastsq(residuals,5,args=(p0,actual), full_output=True) print("actual kd was",kd,"\n") # original value print('\nkd guess', kd_match) # this is the guess for kd # print('cov\n',cov) # inversion of the Hession # print('infodict\n',infodict) # various other outputs print('mesg\n',mesg) # a string with status print('ier\n',ier) # status flag leastsq Example: Graphs m8_leastsq.py (continued): plt.plot(p0,actual) # plot the noisy signal plt.plot(p0,func(kd_match,p0)) # along with the estimate plt.xlabel('p0') plt.ylabel('func(p0)') plt.title('noisy signal and estimate') plt.show() plt.plot(p0,func(kd_match,p0),label='estimate') # plot the estimate plt.plot(p0,clean,label='original') # with clean signal plt.xlabel('p0') plt.ylabel('func(p0)') plt.legend() plt.title('clean original and estimate') plt.show() leastsq Example: Graphs Continued m8_leastsq.py (continued): # difference between predicted value and received value resid = residuals(kd_match,p0,actual) # calculate the residuals plt.plot(p0,resid) # and plot them plt.xlabel('p0') plt.ylabel('residual') plt.title('predicted vs received residuals') plt.show() # difference between the predicted value and the clean value resid = residuals(kd_match,p0,clean) # calculate the residuals plt.plot(p0,resid) # and plot them plt.xlabel('p0') plt.ylabel('residual') plt.title('predicted vs clean residuals') plt.show() leastsq Example: Run Ø python m8_leastsq.py actual kd was 3.0 kd guess [3.25122666] mesg Both actual and predicted relative reductions in the sum of squares are at most 0.000000 ier 1

Use Quizgecko on...
Browser
Browser