• Không có kết quả nào được tìm thấy

Programming for Computations – MATLAB/Octave

Nguyễn Gia Hào

Academic year: 2023

Chia sẻ "Programming for Computations – MATLAB/Octave"


Loading.... (view fulltext now)

Văn bản

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.

Fig. 1.1 Plot generated by the script ball_plot.m , showing the vertical position of the ball at a thousand points in time
Fig. 1.1 Plot generated by the script ball_plot.m , showing the vertical position of the ball at a thousand points in time

More Basic Concepts

  • Using Matlab Interactively
  • Arithmetics, Parentheses and Rounding Errors
  • Variables
  • Formatting Text and Numbers
  • Arrays
  • Plotting
  • 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.

Fig. 1.2 Generated plot for the heights of family members from two families
Fig. 1.2 Generated plot for the heights of family members from two families


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 Tests

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.

For Loops

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.

While Loops

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).

Fig. 3.4 Illustration of the rectangle method with evaluating the rectangle height by either the left or right point
Fig. 3.4 Illustration of the rectangle method with evaluating the rectangle height by either the left or right point

Population Growth

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.

Numerical Solution

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.

Fig. 4.2 Illustration of a forward difference approximation to the derivative
Fig. 4.2 Illustration of a forward difference approximation to the derivative

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/.

Fig. 4.4 Evolution of a population computed with time step 0.5 month
Fig. 4.4 Evolution of a population computed with time step 0.5 month

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

Figure 4.7 shows the resulting curve. We see that the population stabilizes around M D 500 individuals
Figure 4.7 shows the resulting curve. We see that the population stabilizes around M D 500 individuals

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.

Fig. 4.9 Natural evolution of a flu at a boarding school
Fig. 4.9 Natural evolution of a flu 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.

Time-Restricted Immunity

Incorporating Vaccination

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.

Fig. 4.11 Including loss of immunity
Fig. 4.11 Including loss of immunity

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/.

Fig. 4.15 Sketch of a one-dimensional, oscillating dynamic system (without friction)
Fig. 4.15 Sketch of a one-dimensional, oscillating dynamic system (without friction)

Numerical Solution

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.

Fig. 4.16 Simulation of an oscillating system
Fig. 4.16 Simulation of an oscillating system

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.

Fig. 4.19 Illustration of a backward difference approximation to the derivative
Fig. 4.19 Illustration of a backward difference approximation to the derivative

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.

Fig. 4.20 Illustration of a centered difference approximation to the derivative
Fig. 4.20 Illustration of a centered difference approximation to the derivative

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.

Fig. 4.24 The last 10 of 40 periods of oscillations by the 4th-order Runge-Kutta method
Fig. 4.24 The last 10 of 40 periods of oscillations by the 4th-order Runge-Kutta method

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.

Fig. 4.28 Effect of linear damping in combination with a sinusoidal external force
Fig. 4.28 Effect of linear damping in combination with a sinusoidal external force

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
  • Vectorization
  • 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.

Fig. 5.1 Unstable simulation of the temperature in a rod
Fig. 5.1 Unstable simulation of the temperature in a rod


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

Newton’s Method

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.

Fig. 6.1 Illustrates the idea of Newton’s method with f .x/ D x 2  9, repeatedly solving for crossing of tangent lines with the x axis
Fig. 6.1 Illustrates the idea of Newton’s method with f .x/ D x 2 9, repeatedly solving for crossing of tangent lines with the x axis

Making a More Efficient and Robust Implementation

The Secant Method

The Bisection Method

Rate of Convergence

Solving Multiple Nonlinear Algebraic Equations

Abstract Notation

Taylor Expansions for Multi-Variable Functions

Newton’s Method



Error messages

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).

Introducing errors

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.

Fig. 3.1 The integral of v.t / interpreted as the area under the graph of v
Fig. 3.1 The integral of v.t / interpreted as the area under the graph of v


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

Game 21

Linear interpolation

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

Rectangle methods

Adaptive integration

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).

Fig. 4.1 Mesh in time with corresponding discrete values (unknowns)
Fig. 4.1 Mesh in time with corresponding discrete values (unknowns)

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).

Explore symmetry

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

Hình ảnh

Fig. 1.1 Plot generated by the script ball_plot.m , showing the vertical position of the ball at a thousand points in time
Fig. 1.2 Generated plot for the heights of family members from two families
Fig. 3.1 The integral of v.t / interpreted as the area under the graph of v
Fig. 3.2 Computing approximately the integral of a function as the sum of the areas of the trape- trape-zoids

Tài liệu tham khảo

Tài liệu liên quan

Decide whether the following statements are true (T) or false (F). ______ According to Professor Davies, the level of literacy in sixteen-century England matched his