Numerical Methods An overall goal with this book is to motivate computer programming as a very powerful tool for doing mathematics. Comparing these two versions of the book provides an excellent demonstration of how similar these languages are.
What Is a Program? And What Is Programming?
With some programming skills, you might be able to write your own little program that can translate one data format to another. Well, you can write down the recipe in those three languages and forward it.
A Matlab Program with Variables
- The Program
- Dissection of the Program
- Why Not Just Use a Pocket Calculator?
- Why You Must Use a Text Editor to Write Programs
- Write and Run Your First Program
You must know the consequences of every instruction in the program and be able to determine the consequences of the instructions. You see more of these comments in the code, and you probably find that they make it easier to understand (or guess) what the code means.
A Matlab Program with a Library Function
This means that where we see atan(y/x), a calculation is performed (tan1.y=x/) and the result "replaces" the textatan(y/x). With the missing semicolon, Matlab will do the calculations and print the result on the screen.
A Matlab Program with Vectorization and Plotting
This is actually just as magical as if we had written justy/x: then the computation of y/x would take place and the result of that division would replace texty/x. It builds on the example above, but is much simpler, both in terms of the math and number of numbers involved.
More Basic Concepts
- Using Matlab Interactively
- Arithmetics, Parentheses and Rounding Errors
- Formatting Text and Numbers
- Error Messages and Warnings
- Input Data
- Symbolic Computations
- Concluding Remarks
A special program (debugger) can be used to help check (and run) various things in the program that you need to fix. A useful test can then be to remove, say, the last half of the program (by inserting the % comment characters) and insert print commands in smart places to see what happens, for example.
Write a program that calculates the volume of a V-cube with sides LD 4 cm and prints the result on the screen. Both V and L should be defined as separate variables in the program. Then have the program calculate the product of these two variables and print the result on the screen as
If the answer to the "if" question is positive (true), we are done and can skip the next if questions. If the condition is , the following statements up to and including the next if else are executed, and the remaining other branches are skipped.
Notice the two return values result1and result2 specified in the function header ie. the first line of the function definition. This legend (sometimes known as the adoc string) must be placed right at the top of the function.
It's a good rule to develop a program with many functions and then in a later optimization phase, when everything is calculated correctly, remove function calls that have been quantized to slow down the code. At the same time, however, it is something that programmers often do, so it is important to develop the right skills in this area.
Note that the programmer introduced a variable (the loop index) named i, initialized it (i = 1) before the loop, and updated it (i = i + 1) in the loop. In those accidental (incorrect) cases, the boolean expression of the while test never evaluates to false and the program cannot escape the loop.
Reading from and Writing to Files
Compared to aforloop, the programmer does not have to specify the number of iterations when coding awhileloop. If you accidentally enter an infinite loop and the program just hangs forever, press Ctrl+c to stop the program.
Have the program make a formatted printout of the array to the screen, both before and after sorting. Letter numberiis is then reached by gene(i), and a substring of indexiup to and includingj, bygene(i:j) is created. a) Write a function freq(letter, text) that returns the frequency of the letter in the string text, that is, the number of occurrences of letter divided by the length of text.
Basic Ideas of Numerical Integration
If we relax the requirement that the integral be exact, and look instead for approximate values, produced by numerical methods, integration becomes a very simple task for any given .x/(!). That is, f .x/ is very rarely exact, and then it does not make sense to calculate the integral with a smaller error than the one already present inf .x/.
The Composite Trapezoidal Rule
The General Formula
We start with (3.5) and approximate each integral on the right-hand side with a single trapezoid. Strictly speaking, the writing of e.g. "the trapezoidal method" implies the use of only a single trapezoid, while "the compound trapezoidal method" is the most correct name when several trapezoids are used.
Now we compute our special problem by calling application() as the only statement in the main program. The application function and its calls are in the trapezoidal_app.m file, which can be run as.
Alternative Flat Special-Purpose Implementation
Integrand3t2et3 is inserted many times in the code, which quickly leads to errors. How much do we need to change in the previous code to calculate the new integral.
The Composite Midpoint Method
The General Formula
Comparing the Trapezoidal and the Midpoint Methods
The different methods differ in the way they construct the evaluation pointsxi and the weightswi.
- Problems with Brief Testing Procedures
- Proper Test Procedures
- Finite Precision of Floating-Point Numbers
- Constructing Unit Tests and Writing Test Functions
In the trapezoidal and midpoint rules, the error is known to depend on n as much as n2 as n. Solving a problem without numerical errors We know that the trapezoidal rule is exact for linear integrands.
Vectorization essentially eliminates this loop in Matlab (i.e. the overx loop and application to each x value are instead executed in a library of fast, compiled code). Note the need for the vectorized operator. *in the function expression sincev(x) is called with array argumentsx.
Measuring Computational Speed
Double and Triple Integrals
The Midpoint Rule for a Double Integral
2 Cj hy/ : (3.25) Direct derivation Formula (3.25) can also be derived directly in the two-dimensional case using the idea of the method of the mean. The rule of the middle is exact for linear functions, no matter how many subintervals we use.
The Midpoint Rule for a Triple Integral
The derivation of the double integral formula and the implementations follow exactly the same ideas as we explained with the midpoint method, but more terms need to be written into the formulas. Implementation Let's follow the ideas for implementations of the central rule for the double integral.
Monte Carlo Integration for Complex-Shaped Domains
The correct answer is 3, but Monte Carlo integration is unfortunately never exact, so it is impossible to predict the result of the algorithm. Mathematically, it is known that the standard deviation of the Monte Carlo estimate of the integral converges as n1=2, where n is the number of samples.
1as !0, indicating a possible error size problem. pxdx shows that the rate of convergence is actually restored to 2. A remarkable property of the trapezoidal rule is that it is exact for integrals of R sinnt dt (when the subintervals are of equal size).
Derivation of the Model
We are not concerned with the spatial distribution of animals, only with the number of them in a certain space, where there is no exchange of specimens with other spaces. We also present Dbd, which is the net population growth rate per unit of time.
Note that this is an approximation, because the differential equation is originally valid at all real values. Such an algorithm is called a numerical scheme for the differential equation and is often written compactly as.
Programming the Forward Euler Scheme; the Special Case 94
The good thing about Forward Euler's method is that it provides an understanding of what a differential equation is and a geometric picture of how to construct a solution. We know that the line must pass through the solution attn, i.e. point.tn; un/.
Programming the Forward Euler Scheme; the General Case 97
You are now encouraged to do exercise 4.1 to become more familiar with the geometric interpretation of the Forward Euler method. The reader is strongly encouraged to repeat the steps in the derivation of the Forward Euler scheme and establish that we get
Verification: Exact Linear Solution of the Discrete Equations 101
Currently, world population projections point to growth to 9.6 billion before declining. To derive such a model, we can mainly use intuition, so no specific background knowledge of diseases is required.
Spreading of a Flu
The expected number of individuals in category S who catch the virus and become infected in the time interval t is then ptSI. Since there is no loss in the R category (people are either healed and immune or dead), we are done modeling this category.
A Forward Euler Method for the Differential Equation
This differential equation model (and also its discrete counterpart above) is known as the anSIR model. The input data to the differential equation model consists of the parameters ˇ and as well as the initial conditions S.0/ D S0, I.0/ D I0, and R.0/DR0.
Programming the Numerical Method; the Special Case
At another school where the disease had already spread, it was observed that at the beginning of a day there were 40 susceptible and 8 infected, while 24 hours later the numbers were 30 and 18 respectively. We can experiment with and see if we get an outbreak of the disease or not.
Outbreak or Not
This program was written to investigate the spread of a flu at said boarding school, and the rationale for the specific choicesˇand reads as follows. We started out by modeling a very specific case, namely the spread of a flu among students and staff at a boarding school.
Abstract Problem and Notation
We try to incorporate this generalization into the model so that the model has a much wider scope than what we set out to do in the beginning. This is the very power of mathematical modeling: by solving one specific case, we have often developed more general tools that can easily be used to solve seemingly diverse problems.
Programming the Numerical Method; the General Case
Recall that the returned fromode_FE contains all components (S, I, R) in the solution vector at all time points. We can check that this relation holds by comparing SnCinCRn with the sum of the initial conditions.
Furthermore, we assume that in a time interval a fraction of the S category is subject to a successful vaccination. The program must store V .t /in an additional arrayV and the plot command must be expanded with more arguments to plotVversustas well.
Discontinuous Coefficients: a Vaccination Campaign
Oscillating One-Dimensional Systems
Derivation of a Simple Model
At x D0 the spring is not stretched, so the force is zero, and x D0 is therefore the equilibrium position of the body. Equation (4.42) is a second-order differential equation, so we need two initial conditions, one for the position x.0/ and one for the velocity x0.0/.
Programming the Numerical Method; the Special Case
Simulating in three periods the cosine function, T D 3P, and choosing such that there are 20 intervals per period, gives DP =20 and a total of Nt DT = tinterval. Figure 4.16 shows a comparison between the numerical solution and the exact solution of the differential equation.
A Magic Fix of the Numerical Method
The standard way to express this scheme in physics is to change the order of the equations. That is, first the velocity and then the position are updated using the most recently calculated velocity. 4.50) in terms of accuracy, then the order of the original differential equations.
The 2nd-Order Runge-Kutta Method (or Heun’s Method) . 122
The solution of the ODE system is returned as a two-dimensional array, where the first column (sol[:,0]) stores u and the second (sol[:,1]) stores v. It just means that we redefine the name inside the function to mean the instantaneous solution to the first component of the ODE system.
The 4th-Order Runge-Kutta Method
Implementation The phases of the 4th order Runge-Kutta method can be easily implemented as a modification of theosc_Heun.pycode. Note that the 4th-order Runge-Kutta method is completely explicit, so there is never any need to solve linear or non-linear algebraic equations, no matter how it looks.
Illustration of Linear Damping
However, the solution of the dimensionless problem is more general: if we have a solutionu.N tNIˇ/, since then we can find the physical solution of a series of problems.
Illustration of Linear Damping with Sinusoidal Excitation . 137
Due to the contact between the body and the plane, a frictional force f .u0/ also acts on the body. To check that the signs in the definition of off are correct, remember that the actual physical force is f and it is positive (ie f < 0) when acting against a body moving with velocity u0 < 0.
A Finite Difference Method; Undamped, Linear Case
It turns out that this method is mathematically equivalent to the Euler-Cromer scheme. Due to the equivalence of (4.76) with the Euler-Cromer scheme, the numerical results will have the same good properties as constant amplitude.
A Finite Difference Method; Linear Damping
In fact, the Euler-Cromer scheme evaluates a nonlinear damping term as f .vn/when computingvnC1, and this is equivalent to using the backward difference above. Hence, the convenience of the Euler-Cromer scheme for nonlinear damping comes at a cost of increasing the overall accuracy of the scheme from second to first order int.
Assume that the initial condition onu0 is nonzero in the finite difference method of Section 4.3.12:u0.0/DV0. The image below shows snapshots from four different times of temperature evolution.
Finite Difference Methods
- Reduction of a PDE to a System of ODEs
- Construction of a Test Problem with Known Discrete
- Implementation: Forward Euler Method
- Application: Heat Conduction in a Rod
- Using Odespy to Solve the System of ODEs
- Implicit Methods
In other words, we have found a model that is independent of the length of the rod and the material it is made of. Figure 5.4 shows a comparison of the length of all the time steps for two values of the tolerance.
You can then compare the number of time steps with that required by other methods. a) The Crank-Nicolson method for ODEs is very popular when combined with diffusion equations. To avoid oscillations, you should have at most twice the stability limit of the Forward Euler method.
Brute Force Methods
Brute Force Root Finding
A brute force algorithm is to run through all points on the curve and check if one point is below the xaxis and if the next point is above the xaxis, or vice versa. The function must takef and an associated intervalŒa; bas input, as well as a number of points (n), and returns a list of all the roots inŒa; b.
Brute Force Optimization
Model Problem for Algebraic Equations
Deriving and Implementing Newton’s Method
In large industrial applications, where Newton's method solves millions of equations simultaneously, one cannot afford to store all the intermediate approximations in memory, so it is important to understand that the algorithm in Newton's method no longer needs to calculate xnwhen xnC1. This speed of the search for the solution is the primary strength of Newton's method compared to other methods.
Making a More Efficient and Robust Implementation
The Secant Method
The Bisection Method
Rate of Convergence
Solving Multiple Nonlinear Algebraic Equations
Taylor Expansions for Multi-Variable Functions
Matlab has a whole range of ready-made code toolboxes dedicated to specific fields in science and engineering. The present introductory book provides only a small part of all the functionality that Matlab has to offer.
Volume of a cube
Area and circumference of a circle
Volumes of three cubes
Average of integers
Interactive computing of volume and area
Update variable at command prompt
Formatted print to screen
Matlab documentation and random numbers
Note that the first function in a file must have the same name as the name of the file (except the extension .m). By using the reserved word global, a variable can also be known outside the function in which it is defined (without passing it as a parameter).
Compare integers a and b
Functions for circumference and area of a circle
Function for area of a rectangle
Area of a polygon
Write a function polyarea(x, y) that takes two coordinate arrays with the vertices as arguments and returns the area. Test the function on a triangle, a square and a pentagon, where you can calculate the area by alternative methods for comparison.
Average of integers
While loop with errors
Area of rectangle versus circle
Find crossing points of two graphs
Sort array with numbers
The leading order terminates the series for the error, i.e., the error to the smallest power is a good approximation of the error. When we have the solution.N x;N t /, the solution with N dimension Kelvin, reflecting the true temperature in our environment, is given by.
Compute combinations of sets
Write instructions that generate a deck, i.e. all combinations CA,C2,C3, and so on, up to SK. B). A vehicle registration number is on Form DE562, where the letters range from A to Z and the numbers from 0 to 9. Write statements that calculate all possible registration numbers and keep them in a list.
Frequency of random numbers
Test straight line requirement
Fit straight line to data
Fit sines to straight line
Thetrialfunction can execute a loop where the user is prompted for thebn .. values at each pass of the loop and the corresponding graph is displayed. Use this to find and print the smallest error and the corresponding values of b1, b2, and b3.
Count occurrences of a string in a string
As we show below, these tolerances depend on the size of the numbers in the calculations. In addition, many readers of the code will also say that the algorithm looks clearer than in the loop-based implementation.
Hand calculations for the trapezoidal method
Hand calculations for the midpoint method
Compute a simple integral
Hand-calculations with sine integrals
Make test functions for the midpoint method
Explore rounding errors with large numbers
The goal of this exercise is to make a filetest_ode_FE.m that uses the ode_FE function in the fileode_FE.mand automatically verifies the implementation of node_FE. a) The solution calculated by hand in Exercise 4.1 can be used as a reference solution. Figure 5.3 shows four snapshots of the scaled (dimensionless) solutionN.x;N t /.N The power of scale is to reduce the number of physical parameters in a problem, and in the present case we have found one single problem that is independent of the material (ˇ) and the geometry (L).
Integrating x raised to x
Integrate products of sine functions
Revisit fit of sines to a function
Minimization of E with respect to b1; : : : ; bN will give us abest approximation, in the sense that we adaptb1; : : : ; bN. such that SN deviates from as little as possible. Use this property to create a function test_integrate_coeff to verify the implementation of integrate_coeffs. e) Implement the choice f .t / D 1t as a Matlab function f(t) and call integrate_coeffs(f, 3, 100) to see what the optimal choice of b1; b2; b3is. f) Make a function plot_approx(f, N, M, filename) where you plotf(t) together with the best approximationSN calculated as above, using M intervals for numerical integration.
Derive the trapezoidal rule for a double integral
Compute the area of a triangle by Monte Carlo integration
The expression on the left is actually the definition of the derivative N0.t /, so we have The probability that people meet in pairs at time T is (using the empirical frequency definition of probability) equal tom=n, i.e. the number of successes divided by the number of possible outcomes.
Geometric construction of the Forward Euler method
The disadvantage of the backward difference compared to the centered difference (4.80) is that it reduces the order of accuracy in the overall scheme from t2 tot. Using the same trick in the finite difference scheme of the second order differential equation, i.e. using the backward difference inf .u0/ makes this scheme as convenient and accurate as the Euler-Cromer scheme in the general nonlinear casemu00Cf .u0/Cs.u/DF.
Make test functions for the Forward Euler method
Implement and evaluate Heun’s method
Find an appropriate time step; logistic model
Find an appropriate time step; SIR model
Model an adaptive vaccination campaign
Make a SIRV model with time-limited effect of vaccination
Refactor a flat program
Simulate oscillations by a general ODE solver
Equip this file with a test function that reads a file with correct values and compares them to those calculated by the theode_FE function. To find the correct values, modify programsc_FE.m to dump thearray into the file, runosc_FE.m, and let the test function read the reference results from that file.
Compute the energy in oscillations
Use a Backward Euler scheme for population growth
Use a Crank-Nicolson scheme for population growth
Understand finite differences via Taylor series
Write the Taylor series for u.tnCt / (about D tn, as indicated above), then solve the expression in terms of tou0.tn/. Identify, on the right-hand side, the finite difference calculus and infinite series. Write the Taylor series for u.tn/aroundtnC12t and the Taylor series for u.tn Ct / aroundtn C12t.
Use a Backward Euler scheme for oscillations
Subtract the two series, solve in terms of u0.tnC12t /, identify the finite-difference approximation and the error terms on the right-hand side, and write the leading-order error term. Can you use the leading-order error terms in a)–c) to explain the visual observations in the numerical experiment in Exercise 4.12?.
Use Heun’s method for the SIR model
However, the ODE is linear, so a backward Euler scheme leads to a system of two algebraic equations for two unknowns:
Use Odespy to solve a simple ODE
Set up a Backward Euler scheme for oscillations
Implement the method, either yourself from scratch or using Odespy (name isodespy.BackwardEuler). Prove that, in contrast to the forward Euler scheme, the backward Euler scheme leads to significant unphysical damping.
Set up a Forward Euler scheme for nonlinear and damped
Discretize an initial condition
The solution of the equation is not unique unless we also prescribe initial and boundary conditions. In addition, the diffusion equation needs one boundary condition at each point on the boundary@˝ of˝.
Simulate a diffusion equation by hand
We can run it with anything we want, its size just affects the accuracy of the first steps. Rather, you should use more efficient storage formats and algorithms tailored to such formats, but this is beyond the scope of the current text.
Compute temperature variations in the ground
Compare implicit methods
Explore adaptive and implicit methods
Investigate the rule
File name: rod_BE_vs_B2Step.m.. b) The methods Backward Euler, Forward Euler and Crank-Nicolson can have a uniform implementation. For D0 we restore the Forward Euler method, D 1 gives the Backward Euler scheme and D 1=2 corresponds to the Crank-Nicolson method.
Compute the diffusion of a Gaussian peak
The approximation error in the rule is proportional to t, except for D1=2 where it is proportional to tot2. Remarks Although the Crank-Nicolson method, or the rule with D 1=2, is theoretically more accurate than the Backward Euler and Forward Euler schemes, it can exhibit non-physical oscillations as in the present example if the solution is very steep.
Vectorize a function for computing the area of a polygon
Applying Gauss's divergence theorem to the integral on the right-hand side and moving the time derivative outside the integral on the left-hand side leads to. If we interpret the PDE in terms of heat conduction, we can simply explain the result: with the Neumann conditions, no heat can escape from the domain, so the initial heat will just be uniformly distributed, but not escape, so the temperature cannot go to zero (or to the scaled and translated temperature u, to be precise).
Our attention will be limited to Newton's method for such systems of nonlinear algebraic equations. When we solve algebraic equations f .x/D0, we often say that the solution x is a root of the equation.
Solve a two-point boundary value problem
Equations that cannot be reduced to one of those mentioned cannot be solved by general analytical techniques, which means that most algebraic equations that arise in applications cannot be handled with pen and paper. Just move all terms to the left and then the formula to the left of the equal sign isf .x/.
Understand why Newton’s method can fail
See if the secant method fails
Understand why the bisection method cannot fail
Combine the bisection method with Newton’s method
Write a test function for Newton’s method
Solve nonlinear equation for a vibrating beam