Calculating the next element using the previous ones in Matlab. Simple calculations in MATLAB


1.Command Window(Command window).

Mathematical expressions are written on the command line after the " >> " prompt. For example,

To perform the action, press the “Enter” key.

By default, the program writes the result to the special variable ans.

To save the result under a different name, use the variable name and the assignment sign “=”, for example

>> z = 1.25 /3.11

Edit in Command Window You can only use the current command line. In order to edit a previously entered command, you need to place the cursor in the input line and use the “ ” or “ ” keys.

If a command ends with “;”, then the result of its action is not displayed on the command line.

The command window can be closed with the “ ” button, and the “ ” button serves to separate the window from the system interface. You can return to the standard window form using the menu:

Main menu → Desktop → Desktop Layout → Default.

You can clear the command window using the menu:

Main menu → Edit → Clear Command Window.

Main menu of the MatLab system.

Main menu MatLab contains the following items:

File(File) – work with files;

Edit(Edit) – editing;

View(View) – window management;

Web– communication with the developer company via the Internet;

Window(Window) – connection with system windows;

Help(Help) – link to the help system MatLab.

MATLAB system toolbar.

The toolbar buttons have the following purposes:

New file(New) – displays file editor windows;

Open file(Open) – opens file download windows;

Cut(Cut) – cuts the selected fragment and places it on the clipboard;

Copy(Copy) – copies the selected fragment to the clipboard;

Paste(Paste) – transfers the selected fragment from the clipboard to the input line;

Undo(Cancel) – cancels the result of the previous operation;

Redo(Repeat) – repeats the result of the previous operation;

Simulink– creates a Simulink model (extensions MatLab);

Help Window(Help) – opens help windows.

4. Output format of calculation results .



When entering real numbers a dot is used to separate the fractional part!

>> s = 0.3467819

The calculation result is displayed in the format short(short notation for a number), which is defined by default. You can change the format to long(long number notation).

>> format long

0.34678190000000

On the list Numerical Format number formats available

short– a short notation of the number;

long– long number notation;

short e– short notation of a number in floating point format;

long e– long record of a number in floating point format;

short g– the second form of a short notation of a number;

long g– the second form of the long notation of a number;

The format for displaying numeric data can be set in the menu File(file) item Preferences(preferences). Go to tab Command Window(command window). In option Text display(text display), in list Numeric format(numeric format) set short g, in option Numeric display(display numbers) set compact. These output formats correspond to the output of numbers in the universal form of five significant figures and with suppression of space between lines.

Basics of calculations in MatLab.

To perform simple arithmetic operations in MatLab operators are used:

· addition and subtraction +, – ;

· multiplication and division *, / ;

· exponentiation ^ .

Some special variables:

ans – the result of the last operation without an assignment sign;

eps – relative error in floating point calculations;

pi – number;

i or j – imaginary unit;

Inf – infinity;

NaN – undefined value.

Some built-in elementary functionsMatLab:

exp(x) – exponent of the number x;

log(x) – natural logarithm of the number x;

sqrt(x) – square root of number x;

abs(x) – modulus of number x;

sin(x), cos(x), tan(x), cot(x) – sine, cosine, tangent, cotangent of the number x;

asin(x), acos(x), atan(x), acot(x) – arcsine, arccosine, arctangent, arccotangent of number x;

sec(x), csc(x) – secant, cosecant of number x;

round(x) – rounding the number x to the nearest integer;

mod(x,y) – remainder of integer division of x by y, taking into account the sign;

sign(x) – return the sign of number x.

Let's calculate the value of the expression

>> exp(–2.5)*log(11.3)^0.3 – sqrt((sin(2.45*pi)+cos(3.78*pi))/tan(3.3))

If the operator cannot be placed on one line, then it is possible to continue entering it in the next line if you indicate the continuation sign “…” at the end of the first line, for example,

>> exp(–2.5)*log(11.3)^0.3 – ...

sqrt((sin(2.45*pi)+cos(3.78*pi))/tan(3.3))

Functions for working with complex numbers:

Entering a complex number

>> z = 3 + 4i

3.0000 + 4.0000i

The functions abs(z), angle(z) return the modulus and argument of a complex number, where , ;

complex(a,b) constructs a complex number from its real and imaginary parts:

conj(z)returns the complex conjugate number;

imag(z), real(z) returns the imaginary and real part of the complex number z.

6. Vectors.

Entering, adding, subtracting, multiplying by a number.

Vector in MatLab formed using the square brackets operator. In this case, the elements of a column vector are separated by a semicolon “;”, and the elements of a row vector are separated by a space “” or a comma “,”.

Let's introduce a column vector.

>> x =

Let's introduce a row vector .

>> y =

To transpose a vector, use the apostrophe “’”:

>> z = y’

To find the sum and difference of vectors, use the “+” and “–” signs:

>> c = x + z

Multiplication of a vector by a number is carried out both on the right and on the left using the “*” sign.

Vectors can be arguments to built-in functions, e.g.

>> d = sin(c)

To refer to the elements of vectors, parentheses () are used, for example,

>> x_2 = x(2)

The last element of the vector can be selected by typing the command

>> X_end = x(end)

From several vectors you can make one, for example

>> r =

1.3 5.4 6.9 7.1 3.5 8.2

The colon character " : " is used to separate multiple elements from a vector, for example

>> w = r(3:5)

The colon character " : " also allows you to replace elements of a vector, for example,

>> r(3:5)= 0

1.3 5.4 0 0 0 8.2

The “:” symbol can also be used to construct a vector, each element of which differs from the previous one by a constant number, i.e. step, for example

>> h =

1 1.2 1.4 1.6 1.8 2

The step can be negative (in this case, the starting number must be greater than the final number).

A step equal to one can be omitted

>> k =

Basic functions for working with vectors.

  • length(x) – determining the length of vector x;
  • prod(x) – multiplication of all elements of vector x;
  • sum(x) – summation of all elements of vector x;
  • max(x) – finding the maximum element of vector x;
  • min(x) – finding the minimum element of vector x.

If you call the min or max function with two output arguments = min(x),

then the first variable is assigned the value of the minimum (maximum) element, and the second variable is assigned the number of this element.

7 Matrices.

Various ways matrix input.

1. A matrix can be entered as a column vector consisting of two elements, each of which is a row vector and is separated by a semicolon. For example, let's introduce the matrix

>> A =

2. The matrix can be entered line by line by executing the sequence of commands:

>> A = .

The algorithm implemented in the \ operation consists of the following steps:

  • 1. In the trivial case, if A is sparse and diagonal (the section Sparse Matrices is devoted to sparse matrices in MATLAB), the solution is found using the simple formula x k = b k /a kk , where k=1,2,...n.

    • 2. If A is a square, sparse and band matrix, then a band matrix solver is used. There may be two options:

    a. If A is a tridiagonal matrix and b is a single column of real numbers, then the system is solved by Gaussian elimination (its operations are performed only on elements on the diagonals). If during the process of elimination, row permutations need to be done to maintain stability, or if the matrix A is not tridiagonal, then the following point works.
    b. If A is a tape with a density of non-zero elements greater than the parameter bandden, by default equal to 0.5, or the conditions of the previous paragraph are not met, then, depending on the type A and b ( double or single) the following LAPACK library procedures are called:

    A and b are real doubles - procedures DGBTRF, DGBTRS
    A and b complex type double - procedures ZGBTRF, ZGBTRS
    A or b real type single - procedures SGBTRF, SGBTRS
    A or b complex type single - procedures CGBTRF, CGBTRS

    The density of non-zero elements is understood as the ratio of the number of non-zero elements in the tape to the number of all elements in the matrix tape, for example, for a 7 by 7 matrix, the template of which is given below, the number of non-zero elements

    nz = 1 (on diagram No.-2) + 6 (on diagram No.-1) + 7 (on diagram No. 0) + 6 (on diagram No. 1) + 1 (on diagram No. 2) + 1 (on diagram No. 3) = 22,

    and the number of all elements in the tape

    band = 5 (on diagram No.-2) + 6 (on diagram No.-1) + 7 (on diagram No. 0) + 6 (on diagram No. 1) + 5 (on diagram No. 2) + 4 (on diagram No. 3) = 33

    and the density of non-zero elements will be 2/3. Since 2/3 > 0.5, a solver for strip matrices will be used
    Parameter bandden, like other parameters that control sparse matrix algorithms in MATLAB, is set using the spparms function.

    The meaning of this step of the algorithm is that solving a system with a tape matrix does not require actions outside the tape. If the feed is quite full, then there won’t be very many actions with zero elements. On the contrary, if the tape is sufficiently sparse, then greater effect approaches given in further steps of the algorithm can bring.

  • 3. If A is an upper or lower triangular matrix, then the inverse substitution method is used, or, accordingly, the direct substitution method, in which the solution component is found from the last (or first equation) and then the found components are substituted into the following equations to find the next solution components systems of equations.
  • 4. If A - can be reduced to a triangular form by permutations, then the corresponding reduction is performed, and then the system of equations is solved as in step 3.

    • 5. If A is a symmetric matrix (or, in the complex case, Hermitian), and its diagonal contains only positive real elements, then depending on whether A is sparse or not, item a) or item b) holds.

    a. If A is not sparse, then an attempt is made to perform the Cholesky decomposition A = R T R followed by solving systems with triangular matrices R T and R: R T y = b and Rx = y. In this case, the following procedures of the LAPACK library are called:

    A and b are real and double - DLANGE, DPOTRF, DPOTRS, DPOCON
    A and b are complex and double - ZLANGE, ZPOTRF, ZPOTRS, ZPOCON
    A and b are real and of type single - SLANGE, SPOTRF, SPOTRS, SDPOCON
    A and b complex and single type - CLANGE, CPOTRF, CPOTRS, CPOCON

    The Cholesky decomposition will fail if the matrix is ​​not positive definite. But there is no preliminary check for positive definiteness, and it is determined precisely during the decomposition process, since quite often it is enough to perform several steps of the Cholesky decomposition to clarify the fact that the matrix is ​​not positive definite (in the process of calculations it becomes necessary to extract the square root of a negative number). If the Cholesky decomposition is not satisfied, then move on to the next point.

    b. If A is sparse, then symmetrical permutations of rows and columns are first performed symmetric algorithm minimum degree (function symmmd) to reduce the filling of the Cholesky expansion factor, i.e. to reduce the number of new non-zero elements arising during the filling process:

    p = symmmd(A) - vector p contains permutations

    R = chol(A(p, p));

    and two systems with Cholesky multipliers are solved, the first with a transposed multiplier and corresponding permutations in the right-hand side vector

    and the second, with an expansion factor and with entering the solution components into the corresponding positions in the vector x
    x(p, :) = R\y

  • 6. If A is a Hessenberg matrix, then it is reduced to an upper triangular matrix (with appropriate modification of the right-hand side) and then the system is solved by substitutions

    • 7. If A is a square matrix that does not satisfy the conditions of paragraphs 1-6, then depending on whether it is sparse or not, the following actions

    a. If A is not a sparse matrix, then using Gaussian elimination with the choice of a leading element (to ensure stability of the decomposition), an LU decomposition of the matrix A = LU is performed, where

    U - upper triangular matrix
    L is a matrix reduced to triangular by permutations of rows

    and the solution to the system Ax = b is found by sequentially solving systems with triangular matrices Ly = b, Ux = y.
    To perform LU decomposition, the following LAPACK library procedures are called:

    A and b are real and double - DLANGE, DGESV, DGECON
    A and b are complex and double - ZLANGE, ZGESV, ZGECON
    A and b are real and of type single - SLANGE, SGESV, SGECON
    A and b are complex and single type - CLANGE, CGESV, CGECON

    b. If A is a sparse matrix, then the columns are rearranged to reduce the filling of the L and U factors in the process of finding the decomposition. Further, by rearranging rows in the process of LU decomposition, computational stability is ensured, after which systems with triangular matrices are solved again. To perform LU decomposition of a sparse matrix, the procedures of the UMFPACK library are used

    8. If the previous points of the algorithm did not work and, therefore, A is not a square matrix, then everything is determined by its sparsity and one of the points below works:

    a. If A is not a sparse matrix, then the QR decomposition AP = QR is performed, where P is the column permutation matrix, Q is the orthogonal matrix (Q T Q = I), and R is upper triangular. If A is size m by n, then Q is size m by m, and R is size m by n. Next, the solution to the system is found as follows:

    x = P*(R \ (Q" * b))

    since from Ax = b and AP = QR it follows: QRx = bP and Rx = Q T bP.

    b. In the case of a sparse and rectangular matrix A, there is no point in calculating the matrix Q, since it will most likely be filled. Therefore, the QR decomposition algorithm calculates c = Q T b (i.e., the modified right-hand side) and solves the system Rx = c with a triangular matrix to obtain a solution to the original system Ax = b.

In addition to all of the above, the above algorithm evaluates the conditionality of the matrix and displays a warning about a possible error in the solution if the matrix is ​​poorly conditioned (see section The influence of the conditionality of a matrix on the accuracy of solving a system with it). The above algorithm contains all kinds of checks that can take a significant amount of time. Extra time. The next section, The linsolve function for solving systems of linear equations, introduces the linsolve function, which allows you to specify the matrix type, thereby reducing computation time.

linsolve function for solving systems of linear equations

Instead of solving a system (or systems with multiple right-hand sides defined in a matrix) of linear algebraic equations using the backslash sign, you can use the linsolve function, which does not do all the matrix checks inherent in the \ operation algorithm (see previous section).

The linsolve function called in its simplest form with two input arguments and one output argument
x = linsolve(A, b)
solves the system Ax = b in one of the ways, depending on whether the matrix is ​​square or not.

  • If A is a square matrix, then its LU decomposition is first calculated and then two systems with triangular matrices L and U are solved.
  • If A is a rectangular matrix, then its QR decomposition is first calculated and then the system with triangular matrices is solved.

(More details about the actions with the resulting expansion factors are written in the section). Moreover, if matrix A is ill-conditioned or A is a matrix of incomplete rank, then a corresponding warning is displayed, for example:

A = ; b = ; x = linsolve(A,b) Warning: Rank deficient, rank = 1, tol = 4.4686e-015. x = 0 0 2.0000

To suppress output of this message The linsolve function should be called in the command window with two output arguments
= linsolve(A, b)
then the resulting solution will be written in x, and in r - either the rank of the matrix (if A is a rectangular matrix), or the reciprocal of the estimate for its condition number (see section), for example

A = ; b = ; = linsolve(A,b) x = -1.0000e+000 9.9999e-001 3.3307e-006 r = 6.9444e-013

The main advantages of the linsolve function over the \ operation appear when specifying the type of matrix of the system of equations. To do this, before calling the linsolve function, you need to fill out a structure with information about the matrix with the following fields

  • SYM - symmetrical;
  • LT - lower triangular;
  • UT - upper triangular;
  • UHESS - Hessenberg;
  • POSDEF - symmetrical, positive definite;
  • RECT - rectangular;
  • TRANSA - whether it is necessary to solve a system with a matrix transposed to a given one.

Each field can be either true or false. Of course, not all combinations of field values ​​are valid; for example, a matrix cannot be triangular and at the same time symmetric and positive definite. The correct combinations of field values ​​are collected in the following table

LT UT UHESS SYM POSDEF RECT TRANSA
truefalsefalsefalsefalsetrue/falsetrue/false
falsetruefalsefalsefalsetrue/falsetrue/false
falsefalsetruefalsefalsefalsetrue/false
falsefalsefalsetruetruefalsetrue/false
falsefalsefalsefalsefalsetrue/falsetrue/false

If the matrix of the system is positive definite, then this must be taken into account when solving, since for positive definite matrices the solution is based on the Cholesky decomposition, which requires fewer operations than the LU decomposition used when solving systems with general square matrices. This is easy to verify using following example, in which a symmetric positive definite matrix is ​​created (a matrix of random numbers is added with the transpose to it and large numbers are added to the diagonal) and the system of equations with this matrix is ​​first solved as a system with a general matrix (opts.SYM = false and opts.POSDEF = false), and then as with a symmetric and positive definite matrix (opts.SYM = true and opts.POSDEF = true).

% set all fields of the ops structure, except SYM and POSDEF opts.TRANSA = false; opts.UHESS = false; opts.RECT ​​= false; opts.UT = false; opts.LT = false; % create a vector of matrix sizes N = 2.^(8:12); % create empty arrays to record solution times TimeGeneral = ; TimePosSym = ; % in a loop we create matrices and compare solution times for n = N % we create a symmetric and positive definite matrix % and the vector of the right side A = rand(n); A = A + A" + 2*n*eye(n); b = sum(A, 2); % solve the system as a system with a general matrix opts.SYM = false; opts.POSDEF = false; Tstart = cputime; x = linsolve(A,b, opts); Tend = cputime; TimeGeneral = ; % solve the system as a system with a symmetry and position matrix opts.SYM = true; opts.POSDEF = true; Tstart = cputime; x = linsolve( A,b, opts); Tend = cputime; TimePosSym = ; end % display time graphs figure loglog(N, TimeGeneral, N, TimePosSym) legend("TimeGeneral", "TimePosSym")

The resulting computational costs of these methods of solving the system are shown in the graph below

Of course, if the matrix is ​​triangular, then it is very important to indicate this, since solving a system with a triangular matrix is ​​performed in O(n 2) operations, and if a solution algorithm designed for a matrix of a general form is applied to a system with a triangular matrix, then it will take about O(n 3) operations.

The linsolve function does not check whether the matrix satisfies the specified properties, which may result in an error. For example, if, when solving a system with the following matrix and the right-hand side

A = ; b = ; set true for the UT field (and set all others to false) opts.UT = true; opts.TRANSA = false; opts.LT = false; opts.UHESS = false; opts.SYM = false; opts.POSDEF = false; opts.RECT ​​= false; then when solving the system, the linsolve function will consider the matrix as an upper triangular one and select the elements it needs from the upper triangle A x = linsolve(A,b, opts) This will result in a solution x = 1 1 1 corresponding to the matrix A = ;

The influence of matrix conditionality on the accuracy of solving a system with it

In this section, we look at how errors in the elements of the right-hand side vector and in the matrix of a system of linear equations Ax = b can affect the solution of that system. Let us first turn to the case of introducing disturbances into the vector of the right side b. So we solve two systems

Moreover, it is assumed that we solve each of the systems exactly, and the main question that arises is how different the solution x of system (1) will be from the solution of system (2) with a perturbation on the right side δb. This is a rather important question, since the elements on the right side can be obtained as the result of some measurements, respectively containing the error δb k in each component b k. I would like that when measuring the same quantity (each time with its own small errors), the corresponding solutions of systems with slightly different right-hand sides also do not differ very much from each other. Unfortunately, this is not always the case. The reason for this is a characteristic of the matrix called its conditionality. This will be discussed further.

First, you need to introduce a measure of the proximity of vectors and x, i.e. error vector. The measure of the magnitude of a vector is the norm (it can be defined in different ways). For now, let us take as the norm the usual Euclidean norm of a vector (the square root of the sum of the squares of its components), i.e.

Accordingly

where n is the number of unknowns and equations in the system. In addition to the vector norm to estimate the magnitude of the vector, we also need a matrix norm to estimate the magnitude of the matrix. Let’s take the matrix norm, which is defined as follows (it’s called the spectral norm):

those. the spectral matrix norm is equal to the square root of the maximum eigenvalue of the matrix A T A. MATLAB has a norm function for calculating the norms of matrices and vectors, which, in particular, can calculate the above norms. Why we chose these particular norms of vectors and matrices is explained in detail in the section, where a few calculations and definitions are given. This is related to the estimate that we will use for the error in solving the system (the derivation of this inequality is also given in the mentioned section):

Here we denote the condition number of the matrix, which is defined as follows:

From the above inequality it follows that the larger the condition number of the system matrix, the greater the relative error in the solution can be with a small error on the right side.

Let's consider classic example ill-conditioned matrix - the Hilbert matrix - and find out how the error on the right side of the system affects the error in the solution. The Hilbert matrix is ​​defined as follows

To create a Hilbert matrix, MATLAB provides the hilb function, whose input argument specifies the size of the matrix. Let's take a small 6 by 6 matrix:

N = 6; H = hilb(n) H = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.3333 0.2500 0.2000 0.1667 0 .1429 0.1250 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909

X = ones(n, 1); b = H*x b = 2.4500 1.5929 1.2179 0.9956 0.8456 0.7365

We see that there is nothing “suspicious” in either the matrix or the vector on the right side; all the numbers are not very different from each other. Now let's form the perturbed right-hand side b + δb by adding small numbers of the order of 10 -5 to the vector b and solve the system with the perturbed right-hand side to obtain the vector .

Delta_b = 1e-5*(1:n)"; x_tilda = H\(b + delta_b) x_tilda = 0.9978 1.0735 0.4288 2.6632 -1.0160 1.8593

It can be seen that the resulting solution is far from the exact one, where all units should be. Let's calculate the relative error in the solution and on the right side (the norm function by default calculates the Euclidean norm of the vector):

Delta_x = x - x_tilda; LEFT = norm(delta_x)/norm(x) LEFT = 1.1475 RIGHT = norm(delta_b)/norm(b) RIGHT = 2.7231e-005

So, the error in the solution is of the order of unity, although the changes on the right side were of the order of 10 -5. This fits perfectly into the above inequality for error. Indeed, let's calculate the condition number cond(H) using the MATLAB function called cond and by default calculates for the spectral norm of the matrix

C = cond(H) c = 1.4951e+007

LEFT ≈ 1.1475 less

C* RIGHT? 1.4951e+07 * 2.7231e-05 ≈ 407

and the inequality is satisfied (even with some margin).

As the size of the Hilbert matrix increases, the error in the solution will only increase (this is easy to check by setting n = 7, 8, ...). Moreover, for n = 12 a message will be displayed stating that the matrix is ​​ill-conditioned and the solution may be incorrect

Warning: Matrix is ​​close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.409320e-017.

As a measure of conditionality, the value chosen here is RCOND equal to one divided by the estimate of the condition number (the condition number is estimated using a fairly fast algorithm that works much faster than cond, which is described in more detail in the help on the rcond function). If the value of RCOND is small, then the matrix is ​​considered ill-conditioned.

The increase in error in the solution as the size of the Hilbert matrix increases is due to the fact that the condition number of the Hilbert matrix grows very quickly with its size. This is easy to verify using simple loop and semilogy functions (the scale along the ordinate axis is logarithmic):

N = 1:20; C = zeros(1, 20); for n = N H = hilb(n); C(n) = cond(H); end semilogy(N, C) grid on, title("cond(H)"), xlabel("n")

By the way, since the cond function finds the condition number using a numerical method (namely, singular expansion for finding singular numbers), the condition number after n = 12 is no longer calculated correctly; in fact, it should grow further, which can be verified using symbolic calculations in MATLAB and operations with a given number of significant figures

N = 1:20; C = zeros(1, N); digits(60) for n = N H = vpa(sym(hilb(n))); % calculation of the Hilbert matrix accurate to the 60th digit sigma = svd(H); % finding the singular values ​​of the Hilbert matrix % calculating the condition number of the Hilbert matrix C(n) = max(double(sigma))/min(double(sigma)); end semilogy(N, C) grid on, title("cond(H)"), xlabel("n")

Let us now consider three important points regarding this example.

The first of them - solving the system Hx = b (with the right-hand side vector corresponding to the solution from units) using the backslash operator will not give exact units, the error will be already in the tenth digit (although in MATLAB all calculations are performed with double precision by default)

Format long H\b ans = 0.99999999999916 1.00000000002391 0.99999999983793 1.00000000042209 0.99999999953366 1.00000000018389

This happens because when calculating the vector b when multiplying the matrix H by a vector of all ones, we have already included some error in it. In addition, rounding errors in the process of solving the system also played a role and poor conditionality of the matrix (even small size) led to such errors in the solution. Indeed, for small matrices with a small condition number, such an effect will not be observed. Let's take a 6 by 6 matrix of random numbers and calculate its condition number

A = rand(n); cond(A) ans = 57.35245279907571

Then we form the right-hand side corresponding to the exact solution from all units

X = ones(n, 1); b = A*x;

and solve the system, resulting in good accuracy

>> A\b ans = 1.00000000000000 1.000000000000000 1.00000000000000 1.00000000000000 1.000000000000000 1.00000000000000

The second important point is related to the determinant of the matrix. Let's calculate the determinant of the 6 by 6 Hilbert matrix with which we solved the system

>> det(hilb(6)) ans = 5.3673e-018

The determinant turns out to be small enough, from which one can make the incorrect conclusion that the source of the problem of a large error in the solution with a small perturbation of the right side of the system is the smallness of the determinant. In reality, this is not the case; if the determinant of the system is small, this does not mean that the matrix is ​​poorly conditioned and, therefore, there is a large error in the solution. Indeed, let us take the matrix

A = ;

Its determinant is very small, about 10 -121

Det(A) ans = 9.9970e-121

(the word “very” is of course conditional, but it is smaller than the determinant of a 6 by 6 Hilbert matrix, which we had problems with when solving the system). However, the smallness of the determinant will not in any way affect the error in the solution when the right side of the system is perturbed, which is easy to show by forming a system with known solution, introducing disturbances into the right-hand side and solving the system:

X = ; b = A*x; delta_b = 1e-5*; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; LEFT = norm(delta_x)/norm(x) RIGHT = norm(delta_b)/norm(b) LEFT = 2.1272e-005 RIGHT = 2.1179e-005

So, the relative error in the solution is almost equal to the relative error on the right side, since the condition number of the above matrix A is slightly greater than 1, namely:

C = cond(A) c = 1.0303

And finally, consider the third question regarding achieving the equal sign in the inequality for an error in the solution

That is, in other words, let's find out: can there be a situation when we make a small relative perturbation of the right side of the system, say 10 -5, the condition number of the system matrix is ​​equal to 10 10, and in the solution a relative error of 10 5 is obtained. You can go through a lot of examples, trying out different matrices, not achieve equality, and erroneously conclude that this is just an overestimate from above for an error in the solution. However, this is not the case, as the following example convinces us of.

in which the relative perturbation of the right side is 10 -5

RIGHT = norm(delta_b)/norm(b) RIGHT = 1.0000e-005

The relative error in solving the system turns out to be 10 5

X = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; LEFT = norm(delta_x)/norm(x) LEFT = 1.0000e+005

And this happens because
1) in this example, the condition number of matrix A is 10 10 ;

C = cond(A) c = 1.0000e+010

2) the inequality for the error in the solution turns into equality.

If we set the format to long, we see that LEFT, RIGHT and the condition number are not exactly 10 5 , 10 -5 and 10 10 , respectively, but this is due to rounding errors. If solved in exact arithmetic, the equality would be exact, which can be shown analytically (see section).

In previous examples, we assumed that the matrix does not contain errors, and errors can only be on the right side of the system of linear equations. In practice, it may turn out that the elements of the system matrix are specified with errors, i.e. instead of the system

|| Ax || V = || b || V ⇒ || A || M || x || V ≥ || b || V

where || || M is some matrix norm. In order to do this, the matrix norm || || M and vector norm || || V must satisfy the following inequality

|| A || M || x || V ≥ || Ax || V

for any matrices A and vectors x. If this inequality holds, then the matrix norm || || M is called agreed upon with vector norm || || V. It is known, for example, that the spectral norm

(square root of the maximum eigenvalue of the matrix A T) is consistent with the Euclidean vector norm

That is why in this section we carried out all experiments using these norms.

Having received the inequality || A || M || x || V ≥ || b || V further note that from Ax = b it follows that . Since the matrix A is non-singular, it follows that δx = A -1 δb and || δx || V = || A -1 δb || V. Again we use the consistency of norms and obtain the inequality || δx || V ≤ || A -1 || M || δb || V. Further, in the two inequalities obtained

|| A || M || x || V ≥ || b || V and || δx || V ≤ || A -1 || M || δb || V

We divide the smaller value of one of the inequalities by the larger value of the other inequality, and the larger value, respectively, by the smaller one.

and by a simple transformation we finally obtain the required inequality

where cond(A) = || A || M* || A -1 || M.

Note that when deriving the inequality for the error, we used only the fact that the matrix norm is consistent with the vector norm. This is true not only for the spectral norm of a matrix and the Euclidean norm of a vector, but also for other norms. So, for example, the maximum column matrix norm, calculated by the formula

consistent with the first vector norm

and the maximum row matrix norm, calculated by the formula

consistent with the infinite vector norm

The norm function calculates not only the Euclidean norm of a vector and the spectral norm of a matrix, but also the vector and matrix norms listed above. To do this, you need to call it with an additional second argument:

  • q = norm(A, 1) - maximum column norm of matrix A
  • q = norm(A, inf) - maximum row norm of matrix A
  • q = norm(a, 1) - first vector norm of a
  • q = norm(a, inf) - infinite vector norm a

Matrix condition number cond(A) = || A || M* || A -1 || M with respect to various matrix norms can be calculated using the cond function. If cond is called with one input argument (a matrix), then the condition number is calculated relative to the spectral matrix norm. Calling the cond function with an additional argument returns condition numbers relative to the specified matrix norm:

  • c = cond(A, 1) - condition number relative to the maximum column norm of the matrix
  • с = cond(A, inf) - condition number relative to the maximum row norm of the matrix

As an example demonstrating the importance of consistency of the matrix norm with the vector norm in the error inequality, we give an example with matrix A, the right-hand side vector b and the error vector on the right side delta_b.

A = ; b = [ -5.7373057243726833e-001 -1.5400413072907607e-001 -5.3347697688693385e-001 -6.0209490373259589e-001]; delta_b = [-0.71685839091451e-5 0.54786619630709e-5 0.37746931527138e-5 0.20850322383081e-5];

Let us calculate the right and left parts of the estimate using the first vector norm, and the condition number of the matrix in relation to the spectral matrix norm

RIGHT = norm(delta_b,1)/norm(b,1); c = cond(A); x = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tilda; LEFT = norm(delta_x,1)/norm(x,1); format short e disp()

We obtain the following values ​​for the left and right sides of the inequality
9.9484e+004 9.9323e+004

The left side should not exceed the right (as proven above), but it turned out to be larger due to the choice of inconsistent norms.

Let us now examine why the condition number of a matrix cannot be less than one. The condition number of the matrix A is defined as cond(A) = || A || M* || A -1 || M , where || A || M is some matrix norm of A. The matrix norm (i.e. the rule by which a number is associated with each matrix) cannot be arbitrary; it must satisfy the following four axioms:

A1. || A || ≥ 0 for any matrix A and || A || = 0 if and only if A is a zero matrix.

A2. || αA || = | α | * || A || for any matrix A and number α.

A3. || A + B || ≤ || A || + || B || for any matrices A and B

A4. || AB || ≤ || A || * || B || for any matrices A and B.

From the last axiom it is clear that the norm is defined only for square matrices(although in the above formulas for calculating various norms, in principle, there is no such restriction). In addition, from the last axiom it follows that any norm of the identity matrix I is not less than one, indeed

|| I || = || I*I || ≤ || I || 2 ⇒ || I || ≥ 1.

Then, again using the fourth axiom, we find that the condition number of the matrix is ​​always greater than one (true for the condition number of the matrix with respect to an arbitrary matrix norm)

1 ≤ || I || = || AA -1 || ≤ || A || * || A -1 || = cond(A).

Let us complete this section by examining the reason for the appearance of exact equality in the estimate of the error in the solution when the right side is perturbed:

and building corresponding examples in MATLAB. (The discussion below is contained, for example, in the book by J. Forsythe, K. Mohler. Numerical solution of systems of linear algebraic equations. M: “Mir”, 1969.)

An essential role in the reasoning is played by the theorem on the singular decomposition of a matrix, according to which for any real matrix of size n there exist two orthogonal matrices U and V of size n by n (U T U=UU T and V T V = VV T) such that the product D = U T AV is diagonal matrix, and we can choose U and V so that

where μ 1 ≥ μ 2 ≥…≥μ r ≥ μ r+1 =…=μ n =0,

and r is the rank of the matrix A. The numbers μ k are called the spectral numbers of the matrix A. For a non-singular matrix A the following is true:

μ 1 ≥ μ 2 ≥ … ≥μ n ≥ 0.

The next important fact is that multiplication by an orthogonal matrix does not change the Euclidean norm of the vector, i.e. for any vector x with n elements and any orthogonal matrix U of size n by n the equality is true

|| Ux || = || x ||.

Since multiplication by an orthogonal matrix does not change the spectral norm, then

therefore, for the condition number of the matrix the following equality is true:

We have two systems: Ax = b (with exact solution x) and (with exact solution ). Obviously, the error δx is a solution to a system whose right side is a perturbation δb, i.e. systems Aδx = δb. Let D = U T AV be the singular value decomposition of the matrix A, then UDV T = A due to the fact that U and V are orthogonal matrices. Further

Ax = b ⇒ UDV T x = b.

Let us introduce the notation

x" = V T x, b" = U T b,

then the following systems are equivalent

Ax = b ⇔ Dx" = b"

In a completely similar way, consider the system with respect to the error

Aδx = δb ⇒ UDV T δx = δb

We introduce the notation

δx" = V T δx, δb" = U T δb,

for which the following systems turn out to be equivalent

Aδx = δb ⇔ Dδx" = δb"

Those. we have obtained simple equivalent systems with a diagonal matrix, on the diagonal of which are the spectral numbers of the matrix A.

Let us now choose the right-hand sides of these systems in a special way, namely, let

where β > 0, then the corresponding solution to the system Dx" = b" will be

where μ 1 is the maximum singular value of matrix A. We will do the same with the system Dδx" = δ b", namely, let

where γ > 0, then the corresponding solution to the system Dδx" = δb" will be

From the preservation of the norm of a vector when it is multiplied with an orthogonal matrix, it follows that

β/μ 1 = || x" || = || V T x || = || x || and γ/μ n = || δx" || = || V T δx || = ||δx ||.

In exactly the same way we obtain the equalities:

β = || b" || = || U T b || = || b || and γ = || δb" || = || U T δb || = || δb ||.

and since

then we finally get

So, for each matrix A it is possible to construct a vector on the right side of the system and its perturbation such that the error in the solution will be the product of the condition number of the matrix and the error on the right side. In MATLAB, the corresponding construction can be done without using singular value decomposition (although it can be found using the svd function).

First we set n and get two orthogonal matrices U and V by performing a QR decomposition of n by n matrices from random numbers:

N = 4; = qr(rand(n)); = qr(rand(n));

D = diag();

After this, we construct matrix A using the diagonal matrix D and the orthogonal matrices U and V

A = U*D*V";

The condition number of matrix A will coincide with the condition number of matrix D, which in our example is equal to 10 10

Beta = 1; gamma = 1e-5;

and build vectors b" and δb"

B1 = "; db1 = ";

from which we find the vectors b and δb

X = A\b; x_tilda = A\(b+db);

calculate the left and right sides of the inequality

dx = x - x_tilda; RIGHT = norm(db)/norm(b); LEFT = norm(dx)/norm(x);

and bring them out

Format short e disp()

We get equality

Ministry of Education of the Republic of Belarus

Educational institution

· contour3(X, Y, Z) – constructs contour lines for a surface, obtained by stratifying a three-dimensional figure with a series of cutting planes located parallel to the reference plane of the figure.

The example below shows how the commands described can be used to plot the surface Z = sin(X) cos(X).

>> =meshgrid(-3:0.1:3,-3:0.1:3); Z=sin(X).*cos(X);

>> subplot(3,2,1), plot3(X, Y,Z) % Figure 4.3 (a)

>> subplot(3,2,2), mesh(X, Y,Z) % Figure 4.3 (b)

>> subplot(3,2,3), surf(X, Y,Z) % Figure 4.3 (c)

>> subplot(3,2,4), surfc(X, Y,Z) % Figure 4.3 (d)

>> subplot(3,2,5),meshz(X, Y,Z) % Figure 4.3 (e)

>> subplot(3,2,6),contour3(X, Y,Z) % Figure 4.3 (f)


4.3. Formatting graphs

The Matlab system provides the ability to configure and adjust the properties of graphs both using the graphical window interface and by setting the appropriate graphical commands and parameters. In table 4.3.1. Some simple techniques for formatting graphs are given.

Table 4.2 - Formatting graphs

Action

GUI Tools

Switch to editing mode.

Click the button EditPlot in the graphics window toolbar.

Formatting lines and markers of control points of graphs.

In editing mode, double-click on the graph line with the left mouse button. In the window that appears PropertyEditor–Line set all the necessary line parameters (thickness, style, color, etc.).

plot(X, Y, S), plot3(X, Y, Z, S)

(Description of commands is given in clause 4.1 and clause 4.2)

Formatting graph axes.

In edit mode, double-click on the graph axis. In the window that appears PropertyEditor–Axes set all necessary axes parameters.

You can also title a graph and axis labels using the commands InsertTitle, InsertXlabel,InsertYlabel main menu of the graphics window.

axes– controls the properties of the axes.

grid- turns on and off the coordinate grid.

xlabel(S),ylabel(S),zlabel(S)– sets labels near the axes. Here S is a string constant or variable.

title(S)– displays the chart title

Labeling directly onto the chart.

Click the button InsertText,fix the location of the inscription with a mouse click and enter the desired text.

text(X,Y,S)– adds text to a two-dimensional chart, given by string S, so that the beginning of the text is located at the point with coordinates (X, Y).

text(X,Y,Z,S)- adds text to a 3D chart.

Drawing lines and lines with arrows directly onto the chart

Click one of the buttons InsertArraw or InsertLine. Place the mouse pointer at the desired location on the graph and, while holding down the left mouse button, draw a line.

Building a Legend

Insert, and then the command legend.

legend(S1,S2,S3,...)– adds a legend to the current graph with explanations in the form of lines specified in the list of parameters.

Color scale output

In the main menu of the graphic window, select the command Insert, and then the command Colorbar.

colorbar(‘vert'),colorbar(‘horiz')– displays a vertical or horizontal color scale.

Rotating the graph

Click the button Rotate 3D and rotate the graph using the mouse (can also be used for two-dimensional graphs).

rotate3d– sets the rotation of a three-dimensional figure.

with boundary conditions y(t 0 , t end, p) = y, Where t end, t 0 start and end points of intervals. Parameter t(independent variable) does not necessarily mean time, although most often the solution to the DE is sought in the time domain. The DE system in Cauchy form is written similarly to (1.1), but under y in this case, a column vector of dependent variables is implied. Vector p sets the initial conditions.

To solve the DE of the second and higher order they need to be reduced to a first-order DE system.

There are possible differential equations that are not allowed with respect to the derivative:

F(t, y, dy/dt) = 0. (1.2)

Equations (1.2) usually cannot be reduced analytically to form (1.1). However, the numerical solution does not cause any particular difficulties to determine f(y, t) solve (1.2) numerically with respect to the derivative for given y And t.

ODE solvers

Various numerical methods are implemented in MATLAB to solve ODE systems. Their implementations are named solvers ODU.

In this section, the generic name solver means one of the possible numerical methods for solving an ODE: ode45, ode23, ode113, ode15s, ode23s, ode23t, ode23tb, bvp4c or pdepe.

Solvers implement the following methods for solving DE systems:

Ode45 one-step explicit methods Runge-Kutta 4th and 5th orders modified by Dormand and Prince. This is the classic method recommended for initially trying out a solution. In many cases it gives good results, if the system of equations being solved is not rigid.

Ode23 one-step explicit Runge-Kutta methods of 2nd and 4th orders as modified by Bogacki and Champine. With moderate rigidity of the ODE system and low accuracy requirements, this method can provide a gain in solution speed.

Ode113 multi-step Adams-Bashworth-Moulton method of variable order predictor-corrector class. This is an adaptive method that can provide highly accurate solutions.

Ode15s is a multi-step variable order method (from 1 to 5, default 5) using numerical "backward differentiation" formulas. This is an adaptive method and should be used if the ode45 solver does not provide a solution and the remote control system is rigid.

Ode23s is a one-step method using a modified 2nd order Rosenbrock formula. Can provide high speed of calculations with low accuracy of solving a rigid remote control system.

Ode23t implicit trapezoid method with interpolation. This method gives good results when solving problems that describe oscillatory systems with an almost harmonic output signal. For moderately rigid systems, the DE can provide a highly accurate solution.

Ode23tb implicit Runge Kutta method at the beginning of the solution and a method using 2nd order backward differentiation formulas subsequently. Despite the relatively low accuracy, this method may be more effective than ode15s.

Bvp4c serves for the boundary value problem of remote control systems of the form y′ = f(t, y), F(y(a), y(b), p) = 0 (full form system of Cauchy equations). The problems it solves are called two-point boundary value problems, since the solution is sought by specifying boundary conditions both at the beginning and at the end of the solution interval.

All solvers can solve systems of explicit equations y′ = F(t, y), and to solve rigid systems of equations it is recommended to use only special solvers ode15s, ode23s, ode23t, ode23tb.

Using ODE Solvers

tspan vector defining the integration interval [ t 0 t final]. To obtain solutions at specific points in time t 0 , t 1 , …, t final(arranged in order of decreasing or increasing) must be used tspan = [t 0 t 1 … t final];

y 0 vector of initial conditions;

Options argument produced by the odeset function (another odeget or bvpget (bvp4c only) function allows you to print the options set by default or by the odeset/bvpset function);

p 1, p 2,... arbitrary parameters passed to the function F;

T, Y decision matrix Y, where each row corresponds to the time returned in the column vector T.

Let's move on to a description of the syntax of functions for solving remote control systems (the name solver means any of the functions presented above).

[T,Y]=solver(@ F,tspan,y 0) integrates a remote control system of the form y′ = F(t, y) on the interval tspan with initial conditions y 0 . @F descriptor of an ODE function (you can also specify a function in the form " F"). Each row in the solutions array Y corresponds to the time value returned in the column vector T.

[T,Y]=solver(@ F,tspan,y 0 ,options) gives a solution similar to the one described above, but with options determined by the values ​​of the options argument created by the odeset function. Commonly used parameters include relative error tolerance RelTol (default 1e3) and vector acceptable values absolute error AbsTol (all components default to 1e6).

[T,Y]=solver(@ F,tspan,y 0 ,options p 1 ,p 2...) gives a solution similar to the one described above, passing additional parameters p 1 , p 2 , ... in m-file F whenever it is called. Use options= if no options are specified.

Solution of first order ODE

PROCEDURE FOR PERFORMANCE OF THE WORK

· title page;

· initial data of the option;

· the solution of the problem;

· results of solving the problem.

Example

Find a solution to the differential equation on the segment for which at(1,7) = 5,3.

Create a user function in the Command Window

g=@(x,y);

In the function syntax @(x,y) x independent variable y dependent variable x-cos( y/pi) right side of the remote control.

The solution process is carried out by accessing the solver (solver) in the Command Window using the following operator:

Ode23(g,,);

The construction of a graph with a grid is carried out by the following operators:

The result is shown in Fig. 1.1

Rice. 1.2.1. Visualization of the numerical solution

EXERCISE

1. Find solutions to first-order differential equations , satisfying the initial conditions y(x 0 ) = y 0 on the interval [ a,b].

2. Construct graphs of the function.

Task options.

Option No. y(x 0 )=y 0 [a,b]
y 0 (1,8)=2,6
y 0 (0,6)=0,8
y 0 (2,1)=2,5
y 0 (0,5)=0,6
y 0 (1,4)=2,2
y 0 (1,7)=5,3
y 0 (1,4)=2,5
y 0 (1,6)=4,6
y 0 (1,8)=2,6
y 0 (1,7)=5,3
y 0 (0,4)=0,8
y 0 (1,2)=1,4

Laboratory work No. 2

Solving ODE systems

GOAL OF THE WORK

To form students’ ideas about the use of remote control systems in various fields; instill the ability to solve the Cauchy problem for remote control systems.

PROCEDURE FOR PERFORMANCE OF THE WORK

1. Study the theoretical part. Complete the tasks corresponding to the number of your option and demonstrate them to the teacher.

2. Complete a laboratory report, which should contain:

· title page;

· initial data of the option;

· the solution of the problem;

· results of solving the problem.

Example

Solve the system

using the ode23() solver.

Solution:

1. Create an m-file of a function for calculating the right-hand sides of a remote control in the editor.

Let the name in the file editor be sisdu.m, then the function can look like this:

function z=sisdu(t,y)

z1=-3*y(2)+cos(t)-exp(t);

z2=4*y(2)-cos(t)+2*exp(t);

>> t0=0;tf=5;y0=[-3/17,4/17];

>> =ode23("sisdu",,y0);

>>plot(t,y)

Rice. 1.3.1. Visualization of the numerical solution obtained using the ode23 function.

1. What does it mean to solve the Cauchy problem for a remote control system?

2. What methods exist for solving remote control systems?

EXERCISE

1. Find the solution to the remote control system

satisfying the initial conditions on the interval ;

2. Construct function graphs.

For example, the solution function for the 8th option is given:

function z=ssisdu(t,y)

% option 8

z1=-a*y(1)+a*y(2);

z2=a*y(1)-(a-m)*y(2)+2*m*y(3);

z3=a*y(2)-(a-m)*y(3)+3*m*y(4);

z4=a*y(3)-3*m*y(4);

>> =ode23("ssisdu",,);

>> plot(t,100*y)

Rice. 1.3.2. Visualization of the numerical solution obtained using the ode23 function.

Task options.

Option No. Tasks
a m
0,1 1,2
0,2 1,5
0,3 1,7
0,4 1,9
0,5
0,6 1,9
0,7 2,3
0,8 2,7
0,9
0,1 1,5
0,2 1,1
0,3

Laboratory work No. 3

1.4 ODE solution n-th order

GOAL OF THE WORK

To form students’ ideas about the application of higher order remote control in various fields; instill the ability to solve the Cauchy problem for higher-order differential equations using application programs; develop skills in checking the results obtained.

PROCEDURE FOR PERFORMANCE OF THE WORK

1. Study the theoretical part. Complete the tasks corresponding to the number of your option and demonstrate them to the teacher.

2. Complete a laboratory report, which should contain:

· title page;

· initial data of the option;

· the solution of the problem;

· results of solving the problem.

Example 1.

Solve second order differential equations given initial conditions .

Solution:

First we bring the remote control to the system:

1. Create an m-file of the function for calculating the right sides of the remote control.

Let the file name be sisdu_3.m, then the function can look like this:

function z=sisdu_3(x,y)

z2=6*x*exp(x)+2*y(2)+y(1);

2. Perform the following steps:

>> x0=0;xf=10;y0=;

>> =ode23("sisdu_3",,y0);

>> plot(x,y(:,1))

Rice. 1.4.1. Visualization of the numerical solution obtained using the ode23 function.

SAMPLE QUESTIONS FOR JOB DEFENSE

1. What does it mean to solve the Cauchy problem for higher-order differential equations?

2. How to bring the remote control m-th order to the remote control system?

EXERCISE

1. Find a solution to the differential equation that satisfies the initial conditions on the interval.

2. Construct function graphs.

Task options.

Option No. Tasks
Equations Initial conditions







Laboratory work No. 4 – 5

Dynamic systems (DS)

GOAL OF THE WORK

Introducing students to the basic concepts of DS, their classification, phase space of DS, kinematic interpretation of the DS system, evolution of DS. Equation of motion of a pendulum. Dynamics of the Van der Pol oscillator.

2. Dynamic system (DS) a mathematical object corresponding to real systems (physical, chemical, biological, etc.), the evolution of which is uniquely determined by the initial state. The DS is determined by a system of equations (differential, difference, integral, etc.) that allow the existence of a unique solution for each initial condition over an infinite time interval.

The state of the DS is described by a set of variables chosen for reasons of naturalness of their interpretation, simplicity of description, symmetry, etc. The set of states of the DS forms a phase space, each state corresponds to a point in it, and the evolution is depicted by (phase) trajectories. To determine the proximity of states, the concept of distance is introduced in the DS phase space. A set of states at a fixed moment in time is characterized by a phase volume.

The description of DS in the sense of specifying the law of evolution also allows for great variety: it is carried out using differential equations, discrete mappings, using graph theory, the theory of Markov chains, etc. The choice of one of the description methods specifies the specific type of mathematical model of the corresponding DS.

Mathematical model of DS is considered given if dynamic variables (coordinates) of the system are introduced that uniquely determine its state, and the law of evolution of the state over time is indicated.

Depending on the degree of approximation, different mathematical models can be assigned to the same system. Study real systems follows the path of studying relevant mathematical models, the improvement and development of which is determined by the analysis of experimental and theoretical results and their comparison. In this regard, by a dynamic system we will understand precisely its mathematical model. By studying the same dynamic system (for example, the movement of a pendulum), depending on the degree to which various factors are taken into account, we will obtain different mathematical models.


2024 gtavrl.ru.