Practical is the best, teach you how to implement advanced applied mathematics problems in MATLAB (1)


You can refer to the boutique column about the MATLAB series

MATLAB takes you from entry to master in 30 days

MATLAB in-depth understanding of advanced tutorial (with source code)

Friends you like can subscribe by themselves, and your support is my motivation for constant updates!

Problems of Higher Applied Mathematics (001)


The Fibonacci sequence is known by the formula


Can be generated, where the initial value is:


Try to write a MATLAB function to generate a certain Fibonacci value, request

①The function format is y=fib(k), the k-th term ak can be obtained when k is given and assigned to the y vector;

②Write appropriate statements to check the input and output variables to ensure that the function can be called correctly;

③Write this function by means of recursive calling.


function y=fib(n) if round(n)==n & n>=1       if n>=3         y=fib(n-1)+fib(n-2);      else, y=1;      end else     error(’n must be positive integer.’)  end

For example, n = 10 can find the corresponding term as

>> fib(10)

ans =


Now we need to compare the speed of recursive implementation and the speed of loop implementation

>> tic, fib(20), toc

ans =


elapsed_time = 62.0490

>>tic, a=[1 1]; for i=3:30, a(i)=a(i-1)+a(i-2); end, a(30), toc

ans =


elapsed_time =


It should be pointed out that the recursive call is slower and much slower than the loop statement, so it is not particularly necessary. There is no need to use the recursive call to solve this problem.


Code generation for recursive functions

To generate code for MATLAB® recursive functions, the code generator uses compile-time recursion or runtime recursion.

You can influence whether the code generator uses compile-time recursion or run-time recursion by modifying your MATLAB code.

See Force Code Generator to Use Run-Time Recursion.

You can disable recursion or disable runtime recursion by modifying configuration parameters.

When you use recursive functions in MATLAB code for code generation, certain restrictions must be observed.

See Recursive function limits for code generation.

Compile time recursion

With compile-time recursion, the code generator creates multiple versions of the recursive function in the generated code.

The input of each version has a value or size customized for that version.

These versions are called function specializations .

By viewing the MATLAB Function report or the generated C code, you can know that the code generator uses compile-time recursion.

The following is an example of compile-time recursion shown in the report.


Runtime recursion

With runtime recursion, the code generator generates recursive functions in the generated code.

By viewing the MATLAB Function report or the generated C code, you can know that the code generator uses runtime recursion.

The following is an example of runtime recursion shown in the report.


No recursion

In the model configuration parameters, set the Compile-time recursion limit for MATLAB functions to 0.

Disable runtime recursion

Some coding standards (such as MISRA®) do not allow recursion. To increase the likelihood of generating MISRA C®-compliant code, disable runtime recursion.

In the model configuration parameters, clear the Enable run-time recursion for MATLAB functions check box.

If your code requires runtime recursion but runtime recursion is disabled, you must rewrite the code so that it can use compile-time recursion or not.

Recursive function limits for code generation

When using recursion in MATLAB code for code generation, please observe the following restrictions:

The top-level function in the MATLAB Function block cannot be a recursive function, but it can call a recursive function.

Before the first recursive call in a run-time recursive function, assign values ​​to all outputs of the function.

Assign values ​​to all elements of the cell array output of the recursive function at runtime.

The input and output of a recursive function at runtime cannot be a class.

For runtime recursion, the Maximum stack size parameter will be ignored.

Problems of Higher Applied Mathematics (002)


According to matrix theory, if a matrix M can be written as M = A + BCBT, and A, B, C are the corresponding order matrices, then the inverse matrix of the M matrix can be obtained by the following algorithm


Try to write a function to invert the matrix M using MATLAB statements according to the above algorithm, and use a small example to test the program, and compare the accuracy with the direct inversion method.


Write this function

function Minv=part_inv(A,B,C) Minv=inv(A)-inv(A)*B*inv(inv(C)+B’*inv(A)*B)*B’*inv(A);

Suppose the matrix is


And it is known that the matrix can be decomposed into


For this example.


>> M=[51 50 36 16; 50 77 60 32; 36 60 87 48; 16 32 48 68];

iM=inv(M);% numerical inverse, direct solution

iM =

0.0553 -0.0389 0.0017 0.0041

-0.0389 0.0555 -0.0210 -0.0021

0.0017 -0.0210 0.0328 -0.0137

0.0041 -0.0021 -0.0137 0.0244

>> A=diag([1 2 3 4]); B=hankel([1 2 3 4]); C=diag([4 3 2 1]);

iM1=part_inv(A,B,C)% Solving method of block matrix

iM1 =

0.0553 -0.0389 0.0017 0.0041

-0.0389 0.0555 -0.0210 -0.0021

0.0017 -0.0210 0.0328 -0.0137

0.0041 -0.0021 -0.0137 0.0244

At first glance, the results seem to be exactly the same. In fact, the numerical algorithms are different.

Here we use the analytical method to get the inverse of the matrix, and then use the following statement to compare the accuracy of the two results

>> M1=sym(M); iM0=inv(M1)

iM0 =


>> norm(double(iM0)-iM)% The error norm of the direct solution

ans =


>> norm(double(iM0)-iM1)% error norm of indirect solution

ans =


It can be seen that the inverse matrix error obtained by the indirect method is larger, because the inv() function is used many times in the newly written function here, which results in a large transmission error.

It can be concluded from this:

If there is a direct solution to a problem, try not to use indirect methods to increase the transmission error.

Problems of Higher Applied Mathematics (003)


An iterative model is given below


Write the M-function to solve the model. If the initial value of the iteration is x0 = 0, y0 = 0, then please perform 30000 iterations to find a set of x and y vectors, and then light up all the xk and yk coordinates One point (be careful not to connect), and finally draw the desired figure.


The graph drawn in this way is also called the Henon gravitational line graph, which attracts the random points iterated together, and finally obtains a seemingly coherent gravitational line graph.


To solve this problem in a circular form, the Henon gravitational line diagram as shown in the figure can be obtained.

x=0; y=0; for i=1:29999      x(i+1)=1+y(i)-1.4*x(i)^2;      y(i+1)=0.3*x(i); end  plot(x,y,’.’) 

The above algorithm dynamically defines x and y, so it needs to re-dimension every step of the loop, which is time-consuming, so in order to speed up, you can consider pre-defining these two variables, such as x=zeros(1, 30000).



About Henon mapping:

Hénon map is a discrete dynamic system that can produce chaotic phenomena. The iterative expression is:

In the classic Ernon map, the parameter values ​​are respectively taken as a=1.4 and b=0.3. At this time, the system exhibits chaos.

When a and b take other different values, the system can behave as a chaotic phenomenon, an intermittent phenomenon, or converge to a periodic point.

The behavior characteristics of the system under different parameters can be seen through the orbit diagram.

The Ernon map was proposed by the French mathematician Michel Ernon as a simplified model of the Poincaré section of the Lorenz model.

For the classical Ernon map, any initial point either tends to the Ernon's singular attractor or diverges to infinity.

The Ernon attractor has a fractal structure, which is continuous in one direction and a Cantor set in the other direction.

Numerical calculations show that the correlation dimension of the classical Ernon attractor is 1.25±0.02, and the Hausdorff dimension is 1.261±0.003.

Problems of Higher Applied Mathematics (004)


The basic sentence of MATLAB language can obviously draw an equilateral triangle immediately. Try to combine the loop structure and write a small program to draw a series of triangles obtained by rotating the equilateral triangle around its center in the same coordinate system. You can also adjust Rotate the step to observe the effect.


Assuming that the equilateral triangle is rotated by θ degrees counterclockwise, the schematic diagram shown in Figure (a) can be obtained. The three vertices of the triangle are (cos θ, sin θ), (cos(θ + 120◦), sin(θ + 120◦)), (cos(θ + 240◦), sin(θ + 240◦)), the curve can be drawn, as shown in figure (b), try to reduce the step distance, such as choosing ∆θ = 2, 1, 0.1, observe the effect.


 t=[0,120,240,0]*pi/180;% 变换成弧度 xxx=[]; yyy=[];    for i=0:5:360      tt=i*pi/180;     xxx=[xxx; cos(tt+t)];      yyy=[yyy; sin(tt+t)];  end   plot(xxx’,yyy’,’r’),  axis(’square’) 


1. In matlab, the relevant trigonometric functions for calculating angles are sind, cosd, tand, etc.,

In matlab, the relevant trigonometric functions for calculating radians are sin, cos, tan, etc.

2. You can look at the sind function introduction in the help, and enter "help sind" in the command line window.

>> help sin

sin Sine of argument in radians.

sin(X) is the sine of the elements of X.

>> help sind

sind Sine of argument in degrees.

sind(X) is the sine of the elements of X, expressed in degrees.

For integers n, sind(n*180) is exactly zero, whereas sin(n*pi) reflects the accuracy of the floating point value of pi.

3. If the sind function is used to express the angle, enter sind(30), where 30 is the angle.

4. If the angle is expressed by the sin function, enter sin(30/180*pi).

5. You can also use deg2rad to convert the angle to radians. Input sin(deg2rad(30)) to get the same result.

Problems of Higher Applied Mathematics (005)


Select the appropriate step to draw the following graph:


Where t ∈ (−1, 1).


Using the ordinary drawing form, choose equal spacing, and get the curve as shown in Figure a, where x = 0 is rough.

>> t=-1:0.03:1; y=sin(1./t); plot(t,y)

Choose the unequal spacing method, you can get the curve shown in Figure b.

>> t=[-1:0.03: -0.25, -0.248:0.001:0.248, 0.25:0.03:1];





The isometric and unequal-distance strategies here are worth learning, let’s take a closer look!

>> help plot

plot Linear plot.

plot(X,Y) plots vector Y and vector X. If X or Y is a matrix, then the vector will be drawn relative to the rows or columns of the matrix (which row prevails). If X is a scalar and Y is a vector, a disconnected line object will be created and drawn vertically as discrete points at X. plot(Y) plots the columns of Y and their indexes. If Y is complex, plot(Y) is equivalent to plot(real(Y), imag(Y)). In all other usages, the imaginary part is ignored. Various line styles, drawing symbols and colors can be obtained through plot(X,Y,S), where S is a string composed of one element in any or all of the following columns:

b blue

. point


g green

o circle

: dotted

r red

x x-mark

-. dashdot

c cyan

+ plus

- dashed

m magenta

* star

(none) no line

y yellow

s square

k black

d diamond

w white

v triangle (down)

^ triangle (up)

<triangle (left)

> triangle (right)

p pentagram

h hexagram

For example,


plots a cyan dotted line with a plus at each data point;


plots blue diamond at each data point but does not draw any line.


combines the plots defined by the (X,Y,S) triples, where the X's and Y's are vectors or matrices and the S's are strings.

For example, plot(X,Y,'y-',X,Y,'go') plots the data twice, with a solid yellow line interpolating green circles at the data points.

The plot command, if no color is specified, makes automatic use of the colors specified by the axes ColorOrder property.

By default, plot cycles through the colors in the ColorOrder property.

For monochrome systems, plot cycles over the axes LineStyleOrder property.

Note that RGB colors in the ColorOrder property may differ from similarly-named colors in the (X,Y,S) triples.

For example, the second axes ColorOrder property is medium green with RGB [0.5 0], while plot(X,Y,'g') plots a green line with RGB [0 1 0].

If you do not specify a marker type, plot uses no marker.

If you do not specify a line style, plot uses a solid line.

plot(AX,...) plots into the axes with handle AX.

plot returns a column vector of handles to lineseries objects, one handle per plotted line.

The X,Y pairs, or X,Y,S triples, can be followed by parameter/value pairs to specify additional properties of the lines.

For example, plot(X,Y,'LineWidth',2,'Color',[.6 0 0]) will create a plot with a dark red line width of 2 points.


x = -pi:pi/10:pi;

y = tan(sin(x))-sin(tan(x));





See also plottools, semilogx, semilogy, loglog, plotyy, plot3, grid, title, xlabel, ylabel, axis, axes, hold, legend, subplot, scatter.

Advanced Applied Mathematics (006)


Select the appropriate range of θ to draw the following polar coordinate graphs:



The method of drawing a polar coordinate curve is very simple. You can draw a polar coordinate graph with polar(θ,ρ), as shown in Figure 2-4. Pay attention to the point operation when drawing graphics:

 t=0:0.01:2*pi;  subplot(221), polar(t,1.0013*t.^2),% (a)  subplot(222), t1=0:0.01:4*pi; polar(t1,cos(7*t1/2))% (b)  subplot(223), polar(t,sin(t)./t)% (c) subplot(224), polar(t,1-(cos(7*t)).^3)


>> help polar

polar Polar coordinate plot.

polar(THETA, RHO) makes a plot using polar coordinates of the angle THETA, in radians, versus the radius RHO.

polar(THETA, RHO, S) uses the linestyle specified in string S.

See PLOT for a description of legal linestyles.

polar(AX, ...) plots into AX instead of GCA.

H = polar(...) returns a handle to the plotted object in H.


t = 0: .01: 2*pi;

polar(t, sin(2*t) .* cos(2*t),'--r');

See also polarplot, plot, loglog, semilogx, semilogy.

Problems of Higher Applied Mathematics (007)


Use a graphical way to find the approximate solution of the simultaneous equations formed by the following two equations.



These two equations should be drawn with the implicit equation drawing function ezplot(), and the intersection point is the solution of the equation, as shown in Figure a.

ezplot(’x^2+y^2-*x*y^2’); hold on ezplot(’x^3-x^2=y^2-y’)

A more accurate value can be obtained by local magnification, as shown in Figure b.


The two intersection points can be accurately read from the graph, (0.4012, −0.8916), (1.5894, 0.8185).

Try to substitute these two points into the original equation for verification.


>> help ezplot

ezplot (NOT RECOMMENDED) Easy to use function plotter

ezplot is not recommended.  Use FPLOT or FIMPLICIT instead.

ezplot(FUN) plots the function FUN(X) over the default domain -2*PI <X <2*PI, where FUN(X) is an explicitly defined function of X.

ezplot(FUN2) plots the implicitly defined function FUN2(X,Y) = 0 over the default domain -2*PI <X <2*PI and -2*PI <Y <2*PI.

ezplot(FUN,[A,B]) plots FUN(X) over A <X <B.

ezplot(FUN2,[A,B]) plots FUN2(X,Y) = 0 over A <X <B and A <Y <B.

ezplot(FUN2,[XMIN,XMAX,YMIN,YMAX]) plots FUN2(X,Y) = 0 over XMIN <X <XMAX and YMIN <Y <YMAX.

ezplot(FUNX,FUNY) plots the parametrically defined planar curve FUNX(T) and FUNY(T) over the default domain 0 <T <2*PI.

ezplot(FUNX,FUNY,[TMIN,TMAX]) plots FUNX(T) and FUNY(T) over TMIN <T <TMAX.

ezplot(FUN,[A,B],FIG), ezplot(FUN2,[XMIN,XMAX,YMIN,YMAX],FIG), or ezplot(FUNX,FUNY,[TMIN,TMAX],FIG) plots the function over the specified domain in the figure window FIG.

ezplot(AX,...) plots into AX instead of GCA or FIG.

H = ezplot(...) returns handles to the plotted objects in H.


The easiest way to express a function is via a string:

ezplot('x^2-2*x + 1')

One programming technique is to vectorize the string expression using the array operators .* (TIMES), ./ (RDIVIDE), .\ (LDIVIDE), .^ (POWER).

This makes the algorithm more efficient since it can perform multiple function evaluations at once.

ezplot('x.*y + x.^2-y.^2-1')

You may also use a function handle to an existing function.

Function handles are more powerful and efficient than string expressions.



ezplot plots the variables in string expressions alphabetically.

subplot(1,2,1), ezplot('1./z-log(z) + log(-1+z) + t-1')

To avoid this ambiguity, specify the order with an anonymous function:

subplot(1,2,2), ezplot(@(z,t)1./z-log(z) + log(-1+z) + t-1)

If your function has additional parameters, for example k in myfun:


function z = myfun(x,y,k)

z = x.^k-y.^k-1;


then you may use an anonymous function to specify that parameter:


See also

ezcontour, ezcontourf, ezmesh, ezmeshc, ezplot3, ezpolar, ezsurf, ezsurfc, plot, vectorize, function_handle.

Problems of Higher Applied Mathematics (008)


Please draw the three-dimensional graphs and contour lines of xy and sin(xy) respectively.


(a) Just give the following commands, and the resulting graphs are shown in Figures a and b.

 [x,y]=meshgrid(-1:.1:1); surf(x,y,x.*y),figure; contour(x,y,x.*y,30) 

(b) Just give the following commands, and the resulting graphs are shown in Figures c and d.

[x,y]=meshgrid(-pi:.1:pi); surf(x,y,sin(x.*y)), figure;  contour(x,y,sin(x.*y),30)


>> help meshgrid

meshgrid Cartesian grid in 2-D/3-D space

[X,Y] = meshgrid(xgv,ygv) replicates the grid vectors xgv and ygv to produce the coordinates of a rectangular grid (X, Y).

The grid vector xgv is replicated numel(ygv) times to form the columns of X. The grid vector ygv is replicated numel(xgv) times to form the rows of Y.

[X,Y,Z] = meshgrid(xgv,ygv,zgv) replicates the grid vectors xgv, ygv, zgv to produce the coordinates of a 3D rectangular grid (X, Y, Z). The grid vectors xgv,ygv,zgv form the columns of X, rows of Y, and pages of Z respectively.

(X,Y,Z) are of size numel(ygv)-by-numel(xgv)-by(numel(zgv).

[X,Y] = meshgrid(gv) is equivalent to [X,Y] = meshgrid(gv,gv).

[X,Y,Z] = meshgrid(gv) is equivalent to [X,Y,Z] = meshgrid(gv,gv,gv).

The coordinate arrays are typically used for the evaluation of functions of two or three variables and for surface and volumetric plots.

meshgrid and NDGRID are similar, though meshgrid is restricted to 2-D and 3-D while NDGRID supports 1-D to ND. In 2-D and 3-D the coordinates output by each function are the same, the difference is the shape of the output arrays. For grid vectors xgv, ygv and zgv of length M, N and P respectively, NDGRID(xgv, ygv) will output arrays of size M-by-N while meshgrid(xgv, ygv) outputs arrays of size N -by-M.

Similarly, NDGRID(xgv, ygv, zgv) will output arrays of size M-by-N-by-P ​​while meshgrid(xgv, ygv, zgv) outputs arrays of size N-by-M-by-P.


Evaluate the function x*exp(-x^2-y^2) over the range -2 <x <2, -4 <y <4,

[X,Y] = meshgrid(-2:.2:2, -4:.4:4);

Z = X .* exp(-X.^2-Y.^2);


Class support for inputs xgv,ygv,zgv:

float: double, single

integer: uint8, int8, uint16, int16, uint32, int32, uint64, int64

See also surf, slice, ndgrid.

>> help surf

surf 3-D colored surface.

surf(X,Y,Z,C) plots the colored parametric surface defined by four matrix arguments. The view point is specified by VIEW. The axis labels are determined by the range of X, Y and Z, or by the current setting of AXIS. The color scaling is determined by the range of C, or by the current setting of CAXIS. The scaled color values ​​are used as indices into the current COLORMAP. The shading model is set by SHADING.

surf(X,Y,Z) uses C = Z, so color is proportional to surface height.

surf(x,y,Z) and surf(x,y,Z,C), with two vector arguments replacing the first two matrix arguments, must have length(x) = n and length(y) = m where [m, n] = size(Z). In this case, the vertices of the surface patches are the triples (x(j), y(i), Z(i,j)).

Note that x corresponds to the columns of Z and y corresponds to the rows.

surf(Z) and surf(Z,C) use x = 1:n and y = 1:m.

In this case, the height, Z, is a single-valued function, defined over a geometrically rectangular grid.

surf(...,'PropertyName',PropertyValue,...) sets the value of the specified surface property. Multiple property values ​​can be set with a single statement.

surf(AX,...) plots into AX instead of GCA.

surf returns a handle to a surface plot object.

AXIS, CAXIS, COLORMAP, HOLD, SHADING and VIEW set figure, axes, and surface properties which affect the display of the surface.

See also surfc, surfl, mesh, shading.

>> help contour

contour Contour plot.

contour(Z) draws a contour plot of matrix Z in the xy plane, with the x-coordinates of the vertices corresponding to column indices of Z and the y-coordinates corresponding to row indices of Z. The contour levels are chosen automatically.

contour(X,Y,Z) draws a contour plot of Z using vertices from the mesh defined by X and Y. X and Y can be vectors or matrices.

contour(Z,N) and contour(X,Y,Z,N) draw N contour lines, choosing the levels automatically.

contour(Z,V) and contour(X,Y,Z,V) draw a contour line for each level specified in vector V. Use contour(Z,[vv]) or contour(X,Y,Z,[vv] ) to draw contours for the single level v.

contour(AX, ...) plots into the axes AX.

[C,H] = contour(...) returns contour matrix C and a handle, H, to a contour object. These can be used as inputs to CLABEL. The structure of a contour matrix is ​​described in the help for CONTOURC.

contour(..., LineSpec) draws the contours using the line type and color specified by LineSpec (ignoring marker symbols).

To specify additional contour properties, you can follow the arguments in any of the syntaxes described above with name-value pairs.


[c,h] = contour(peaks);


See also contour3, contourf, clabel.