Many equations in science have no formula for their solution. In this course the computer solves them. You find a bracket with a change of sign, write the bisection method, and then Newton's method, which doubles the number of correct digits at every step, and learn when each of them fails. You work with polynomials in NumPy, with polyval, roots and polyfit, fit a straight line to measurements and judge the fit, and solve systems of linear equations with matrices, including a circuit with three unknowns. A project uses brentq and curve_fit from SciPy to fit the decay of a radioactive sample, and finds its half-life with an uncertainty. The practice runs your code and tells you which results are wrong.
5 lessons, about 3 h 25 min5 simulations10 problems10 practice questionsPython in every lesson
Many equations have no formula for their solution, such as cos x equal to x, and a computer finds the root, the value where f of x is 0, by calculation. A continuous function that has opposite signs at a and b has a root between them. Bisection takes the middle point, keeps the half in which the sign changes, and repeats, so that the interval holds the root and its width halves at every step. From a width of 0.5, about 19 steps give a width below one millionth. It works for square roots and for the time when a ball lands. It needs a change of sign, so it fails when f has the same sign at both ends, and for a root where the curve only touches the axis, such as x squared at 0. np.sign and np.where can scan an array of values for the places where the sign changes.
Newton's method starts from a guess x, draws the tangent line to the graph at that point, and takes the place where the tangent meets the x axis as the next guess, x minus f of x over the slope f prime of x. Near a simple root the error is squared at every step, so the number of correct digits doubles, and the root of x squared minus 2 is found to full accuracy in 5 steps, where bisection needs about 39 for the same accuracy on a cubic. The method can fail: a start can send it to a cycle, as for x cubed minus 2x plus 2 from 0, and a zero slope makes a division by zero. A safe plan is to find a bracket first, or to use a start that is near the root.
NumPy stores a polynomial as an array of coefficients from the highest power down. np.polyval gives its values at many points, np.roots gives its roots, which may be complex numbers when there is no real root, np.poly builds the coefficients from the roots, and np.polyder gives the derivative. np.poly1d wraps the coefficients in an object that can be called, multiplied and differentiated. np.polyfit with a degree finds the polynomial that fits measured data best by the method of least squares, and with degree 1 it gives the slope and the intercept of a straight line. The residuals are the differences between the data and the fitted values, and R squared tells how much of the scatter of the data the fit explains. A fit is only as good as the model behind it.
A system of linear equations is written with a matrix A of the coefficients, a vector x of the unknowns and a vector b of the right sides, as A times x equal to b. In NumPy, @ is the matrix product, * works place by place, and np.linalg.solve of A and b solves the system. np.allclose of A at x and b checks the answer. Two equations are two lines: they cross at one point when the determinant is not 0, are parallel with no solution, or are the same line with many solutions, and in the last two cases the determinant is 0 and np.linalg.solve raises a LinAlgError. The method works for any number of unknowns, for example the currents of a circuit found with the laws of Kirchhoff.
SciPy has ready tools for the methods of this chapter. scipy.optimize.brentq solves an equation f of x equal to 0 on a bracket with a change of sign, as fast as the best methods and as safe as bisection. scipy.optimize.curve_fit finds the numbers of any model function that fit measured data best, and returns the best values and a covariance matrix whose diagonal gives their uncertainties. With a sigma of the square root of the counts, a fit of the exponential decay n0 times e to the minus lambda t to the counts of a radioactive sample gives the decay constant, and the half-life is ln 2 over lambda, about 20 seconds in the example. The same half-life comes from a straight line through the logarithms of the counts, and brentq finds the time when the counts reach 100.
Log in to join the discussion.
Loading…