Calculs et approximation de données dans MATLAB

Résolution de systèmes d'équations algébriques linéaires et de problèmes spectraux

La résolution de systèmes d'équations algébriques linéaires est l'un des principaux problèmes de calcul, car un grand nombre de problèmes découlant d'une grande variété de domaines appliqués y sont réduits. Les méthodes numériques utilisées pour la solution approximative de problèmes mécaniques, pour le calcul du débit de liquides et de gaz, d'autres processus physiques en milieu continu, les méthodes de calcul de circuits électriques, les méthodes d'approximation des données conduisent à la nécessité de résoudre des systèmes d'équations algébriques linéaires (et il ne s'agit pas d'une liste complète des sources d'occurrence des systèmes algébriques linéaires). Des problèmes spectraux surviennent, par exemple, dans l'analyse de fréquence des structures, dans l'étude de la stabilité des processus dans les circuits électriques et les systèmes de contrôle, ainsi que dans de nombreux autres domaines appliqués.

Présentation des capacités de MATLAB pour résoudre des systèmes d'équations algébriques linéaires

L'environnement MATLAB a été initialement développé pour travailler avec des matrices (MATLAB est l'abréviation de Matrix Laboratory), donc l'arsenal d'outils MATLAB pour résoudre des systèmes d'équations algébriques linéaires est assez riche et comprend :

• résolution de systèmes avec des matrices carrées et rectangulaires ;
• résolution de systèmes par des méthodes directes et itératives (y compris la possibilité de préconditionnement) ;
• décompositions matricielles ;
• stockage de grandes matrices clairsemées sous une forme compacte et algorithmes spéciaux pour résoudre des systèmes avec de telles matrices.

Résolution de systèmes à l'aide de barres obliques inverses

La manière la plus simple de résoudre des systèmes consiste à utiliser le signe barre oblique inverse. Supposons que nous devons résoudre le système

Pour ce faire, remplissez la matrice et le vecteur colonne sur le côté droit (le côté droit doit être une colonne, sinon une erreur sur l'inadéquation des dimensions sera affichée)

et utilisez le signe antislash

X = A\f x = 1 1 1

Au lieu du signe barre oblique inverse, vous pouvez appeler la fonction mldivide

X = mldivide(A, f)

Le résultat sera le même

En calculant l'écart, nous nous assurons que la solution a été trouvée correctement (d'une manière générale, de petites valeurs des composants d'écart n'indiquent pas toujours une solution correctement trouvée)

F - A*x ans = 0 0 0

Il est très important de ne pas confondre A et f, car lors de l'exécution d'une opération

il n'y aura pas d'erreur, mais elle sera affichée

X = 0,2707 0,3439 0,3854

ce qui n'a rien à voir avec la solution du système considéré. Voyons pourquoi cela se produit en écrivant un système avec une « matrice » f et un « côté droit » A. Dans ce cas, il n’y aura qu’une seule inconnue, car f taille 3 par 1 :

S'il y a une matrice sur le côté droit, alors autant de systèmes seront résolus qu'il y a de colonnes de la matrice, et chaque colonne est un vecteur du côté droit du système correspondant (c'est-à-dire que l'opérateur barre oblique inverse vous permet de résoudre plusieurs systèmes à la fois). Ainsi, les systèmes sont résolus

À proprement parler, aucun de ces systèmes n’a de solution. Cependant, si la matrice du système est rectangulaire et que le nombre de ses lignes est supérieur au nombre de colonnes (c'est-à-dire que le nombre d'équations est supérieur au nombre d'inconnues), alors un tel système est appelé systèmes surdéterminés et son approximation la solution est recherchée en minimisant la norme euclidienne du résidu. Donc pour le premier système la solution sera X, qui fournit un minimum à l'expression suivante

(4-7x 1) 2 +(3-11x 1) 2 +(2-12x 1) 2

Il est facile de s’assurer que le minimum est suffisant. Utilisons, par exemple, la Symbolic Math Toolbox

Syms x R = (7*x-4)^2 + (11*x-3)^2 + (12*x-2)^2 dRdx = diff(R) x = résoudre (dRdx) x = 85/314 double(x) ans = 0,2707

Ainsi, le vecteur obtenu en écrivant x = f\A contient la solution de trois systèmes surdéterminés dont le côté droit de chacun est la colonne correspondante de la matrice A.

Bien entendu, la matrice du système surchargé ne contient pas nécessairement une seule colonne. À titre d'exemple, considérons le problème de sélection de paramètres un Et b modèle linéaire selon les données tabulaires données

k	1	2	3	4	5
xk	1.0	1.5	2.0	2.5	3.0
ouais	2.99	2.81	2.89	3.03	3.21

Ce problème se réduit à un système surdéterminé par rapport à un Et b, si nous exigeons l'égalité y(xk)=ouais Pour k=1,2,...,5:

Nous composons la matrice système et le vecteur du côté droit

X = (1:0.5:3)" A = ; f = ; et résolvez-le c = A\f nous obtenons c = 2,9903 2,0100 c'est-à-dire que le graphique des données et des approximations confirme l'exactitude des résultats obtenus figure fun = @(x ) c(1)/x + c(2)*log(x); fplot(fun, ) tenir sur plot(x, y, "ou")

Si un système a plus d'équations que d'inconnues, alors un tel système est appelé système sous-déterminé et peut être résolu dans MATLAB à l'aide du signe barre oblique inverse. Il aura une infinité de solutions et une solution est trouvée qui contient autant de valeurs nulles que possible. Considérons cet exemple

UNE = ; f = ; x = A\f x = 1,3333 1,0000 0 1,6667

En plus du caractère barre oblique inverse, un caractère barre oblique ordinaire est également utilisé. Ils sont reliés par la règle suivante

B/A est équivalent à (A"\B")"
où l'apostrophe signifie transposition. Au lieu de la barre oblique habituelle, vous pouvez également utiliser la fonction mrdivide.

Pour l'instant, nous n'entrerons pas dans les algorithmes de résolution de systèmes intégrés dans l'opération antislash. Cette rubrique y est dédiée. En particulier, la décomposition LU est utilisée pour résoudre des systèmes avec une matrice carrée générale. De plus, la conditionnalité de la matrice est vérifiée. Un exemple classique de matrice mal conditionnée est la matrice de Hilbert, dont les éléments sont déterminés par la formule. Pour créer une matrice de Hilbert d'une taille donnée dans MATLAB, utilisez la fonction hilb. Résolvons, par exemple, un système d'ordre 13 dont le membre de droite est tel que sa solution doit être toutes les unités ( k Le ème élément du côté droit est la somme des éléments k-ième ligne de la matrice). A = hile(13); f = somme(A, 2); x = A\f

En conséquence, nous recevons un message concernant un mauvais conditionnement

Attention : La matrice est proche du singulier ou mal mise à l'échelle. Les résultats peuvent être inexacts. RCOND = 2,339949e-018.

(RCOND est une estimation de l'inverse du numéro de condition) et une solution incorrecte

X = 1,0000 1,0000 1,0004 0,9929 1,0636 0,6552 2,2023 -1,7860 5,3352 -3,4773 3,9431 -0,1145 1,1851

Les dangers générés par un mauvais conditionnement de la matrice sont décrits dans la section. Passons maintenant à un aperçu des décompositions matricielles disponibles dans MATLAB.

Décompositions matricielles Cholesky, LU et QR

MATLAB a des fonctions pour les extensions suivantes :

• Expansion Cholesky - fonction chol ;
• Extension LU - fonction lu ;
• Décompositions QR - fonction qr.

Les extensions aident à résoudre efficacement tout un ensemble de systèmes d'équations algébriques linéaires avec la même matrice et plusieurs vecteurs du membre droit. Considérons l'ensemble de systèmes suivant :

Axe=f(k), k=1, 2,...,N.

Supposons que la matrice UN se présente comme un produit de deux matrices, et des systèmes avec des matrices A=BC et systèmes avec matrices B Et C résolu beaucoup plus rapidement qu'avec une matrice UN. Ensuite, les solutions aux systèmes originaux sont obtenues comme suit. Parce que le Axe=f(k) Et A=BC, Que BCx=f (k) et résoudre séquentiellement Par=f(k) Et Cx=y nous trouvons une solution kème système. Pour les matrices remplies, la décomposition s'effectue comme suit : Ô(n 3) opérations où n- la taille de la matrice (c'est-à-dire en même temps que la résolution du système), et la résolution du système avec chaque facteur d'expansion en Ô(n 2) opérations. Donc la solution N systèmes c avec expansion préalable prend Ô(n 3)+Ô(n 2)N. Résoudre chaque système sans décomposition matricielle nécessiterait Ô(n 3)N opérations arithmétiques.

Décomposition de Cholesky pour une matrice donnée UN consiste à trouver une telle matrice triangulaire supérieure R. avec des éléments diagonaux positifs, qui A=RTR. On sait que si la matrice UN symétrique et défini positif (c'est-à-dire pour tout vecteur X est vrai : x T Ax0, ou, ce qui revient au même, toutes les valeurs propres de la matrice sont positives), alors la décomposition de Cholesky existe et est unique.

Considérons comme exemple la décomposition de Cholesky de la matrice

UNE = ; R = chol(A) R = 2,0000 0,5000 0,5000 0 1,9365 0,3873 0 0 1,8974

Un simple contrôle permet de s'assurer que la décomposition est effectuée correctement (dans la limite des erreurs survenant lors d'opérations sur des nombres réels)

A - R"*R ans = 1,0e-015 * 0 0 0 0 0 0 0 0 0,4441

Si la matrice n'est pas définie positive, cela deviendra clair au cours du processus de décomposition, car il sera nécessaire d'extraire la racine d'un nombre négatif et un avertissement correspondant s'affichera :

UNE = ; R = chol(A)??? Erreur lors de l'utilisation de ==> chol La matrice doit être définie positive.

La fonction chol ne vérifie pas la symétrie de la matrice. Les calculs utilisent des éléments de la matrice originale situés sur la diagonale et au dessus :

UNE = ; R = chol(A) R = 2,0000 0,5000 1,0000 0 1,9365 0,2582 0 0 1,7127 A - R"*R ans = 0 0 0 0 0 0 -1 0 0

Pour une matrice carrée arbitraire, on peut faire L.U.-la décomposition, c'est-à-dire trouver la matrice triangulaire inférieure L et matrice triangulaire supérieure U tel que A = LU. Une déclaration plus rigoureuse est formulée comme suit.

Si Ak- mineur principal d'une matrice carrée UN taille n, compilé à partir du premier k lignes et colonnes (c'est-à-dire A(1:k, 1:k)) et det( Un k)0 pour k=1,2,..., n-1, alors il existe une unique matrice triangulaire inférieure L, dont la diagonale est constituée de uns, et la seule matrice triangulaire supérieure U tel que A = LU.

Lors du calcul L.U. Une expansion LU peut nécessiter une permutation de lignes pour assurer la stabilité numérique du processus d'expansion, afin que la fonction lu puisse renvoyer une matrice L, qui est triangulaire inférieur jusqu'à la permutation des lignes, par exemple :

UNE = ; = lu(A) L = 0 1,0000 0 1,0000 0 0 0,5000 0,5000 1,0000 U = 2 3 1 0 1 1 0 0 4

Travail L.U. Où L- triangulaire jusqu'à permutation des lignes, et sera égal à la matrice UN(d'une manière générale, en tenant compte de l'erreur dans les opérations avec des nombres réels).

A-L*U ans = 0 0 0 0 0 0 0 0 0

Si pour la matrice UN Si la décomposition de Cholesky ou la décomposition LU est effectuée, alors la solution du système avec la matrice UN, comme mentionné ci-dessus, revient à résoudre deux systèmes à matrices triangulaires. Cette solution peut être réalisée en utilisant l'opération de barre oblique inverse, puisque l'algorithme sous-jacent identifie les matrices triangulaires et applique une manière efficace pour résoudre le système avec elles, nécessitant Ô(n 2) opérations arithmétiques. Considérons un exemple de résolution du système

avec préliminaire L.U.-décomposition matricielle.

UNE = ; f = ;

Nous effectuons L.U.-décomposition

Lu(UNE);

et résoudre séquentiellement deux systèmes à matrices triangulaires, d'abord avec L, puis avec U

Y = L\f; x = U\y x = 1 1 1

La solution de deux systèmes peut être écrite en une seule expression

X = U\(L\f)

et le résultat sera le même. Il convient de prêter attention à l'importance d'utiliser des parenthèses pour déterminer l'ordre de solution des systèmes à matrices triangulaires. Expression sans parenthèses

X = U\L\f

conduira à un résultat complètement différent, puisqu'il est exécuté en premier U\L, ce qui équivaut (comme discuté ci-dessus) à la résolution de systèmes avec la matrice U et les colonnes L comme vecteurs de droite. Le résultat est une matrice dont chaque colonne est une solution du système correspondant, puis avec cette matrice et le vecteur de droite f le système est résolu. Ce processus n’a évidemment rien à voir avec la solution du système original.

Considérons la dernière décomposition matricielle - QR-décomposition effectuée par la fonction qr. QR-la décomposition peut être effectuée pour des matrices rectangulaires, précisément si UN matrice de taille m sur n, alors il existe une matrice orthogonale Q taille m sur m(c'est-à-dire tel que Q -1 = Q T) et la matrice R. taille m sur n avec tous les éléments nuls sous la diagonale principale, ce qui A=QR.

QR-l'expansion a une bonne stabilité de calcul et est utilisée en particulier lors de l'approximation de données à l'aide de la méthode des moindres carrés. À titre d'exemple, envisagez de résoudre le système précédent en utilisant la décomposition QR :

Qr(A) Q = -0,7428 0,6000 -0,2971 -0,5571 -0,8000 -0,2228 -0,3714 0,0000 0,9285 R = -5,3852 -5,3852 -5,0138 0 -5,0000 0,4000 0 0 6,6108

Une fois fait QR-décomposition, solution du système Hache=f peut être trouvé comme une solution au système matriciel triangulaire R., parce que parce que le QR=f, Que Rx=QTf:

>> x = R\(Q"*f) x = 1,0000 1,0000 1,0000

Plus de détails sur les décompositions matricielles, les fonctions MATLAB correspondantes chol, lu et qr et leurs arguments sont écrits dans la section Décompositions matricielles.

Algorithme de résolution d'un système d'équations linéaires à l'aide du signe barre oblique inverse

Lors de la résolution d'un système d'équations algébriques linéaires dans MATLAB à l'aide du signe barre oblique inverse
x = UNE\b
fonctionne un algorithme qui, selon le type de matrice, résout le système en utilisant la méthode la plus appropriée, implémentée dans l'une des procédures du package LAPACK ou UMFPACK. Ici b peut aussi être une matrice dont le nombre de lignes coïncide avec le nombre de lignes de la matrice A. Ensuite la matrice x est renvoyée, chaque i-ème colonne contient la solution du système Ax (i) = b (i ) , i=1,2,... ,k, où b (i) est la i-ième colonne de la matrice b = [ b (1) | b (2) ...| b(k)].

L'algorithme implémenté dans l'opération \ comprend les étapes suivantes :

• 1. Dans le cas trivial, si A est clairsemé et diagonal (la section Sparse Matrices est consacrée aux matrices clairsemées dans MATLAB), la solution est trouvée à l'aide de la formule simple x k = b k /a kk , où k=1,2,. ..n.

2. Si A est une matrice carrée, clairsemée et en bandes, alors un solveur de matrice de bandes est utilisé. Il peut y avoir deux options :

un. Si A est une matrice tridiagonale et b est une seule colonne de nombres réels, alors le système est résolu par élimination gaussienne (ses opérations sont effectuées uniquement sur les éléments des diagonales). Si lors du processus d'élimination, pour maintenir la stabilité, il est nécessaire de faire des permutations de lignes, ou si la matrice A n'est pas tridiagonale, alors le point suivant fonctionne.
b. Si A est une bande avec une densité d'éléments non nulle supérieure au paramètre bandé, par défaut égal à 0,5, ou si les conditions du paragraphe précédent ne sont pas remplies, alors, selon le type A et b ( double ou célibataire) les procédures suivantes de la bibliothèque LAPACK sont appelées :

A et b sont de vrais doubles - procédures DGBTRF, DGBTRS
Types complexes A et B doubles - procédures ZGBTRF, ZGBTRS
A ou b type réel unique - procédures SGBTRF, SGBTRS
Type complexe A ou B unique - procédures CGBTRF, CGBTRS

La densité d'éléments non nuls s'entend comme le rapport du nombre d'éléments non nuls dans la bande au nombre de tous les éléments de la bande matricielle, par exemple pour une matrice 7 par 7 dont le modèle est donné ci-dessous, le nombre d'éléments non nuls

nz = 1 (sur schéma n°2) + 6 (sur schéma n°1) + 7 (sur schéma n°0) + 6 (sur schéma n°1) + 1 (sur schéma n°2) + 1 (sur schéma n°3) = 22,

et le nombre de tous les éléments de la bande

bande = 5 (sur schéma n°2) + 6 (sur schéma n°1) + 7 (sur schéma n°0) + 6 (sur schéma n°1) + 5 (sur schéma n°2) + 4 (sur schéma n°3) = 33

et la densité des éléments non nuls sera de 2/3. Puisque 2/3 > 0.5, un solveur pour matrices bandes sera utilisé
Paramètre bandé, comme les autres paramètres qui contrôlent les algorithmes matriciels clairsemés dans MATLAB, est défini à l'aide de la fonction spparms.



La signification de cette étape de l’algorithme est que la résolution d’un système avec une matrice de bande ne nécessite pas d’actions en dehors de la bande. Si le flux est assez plein, il n’y aura pas beaucoup d’actions avec zéro élément. Au contraire, si la bande est assez clairsemée, les approches données dans les étapes ultérieures de l'algorithme peuvent alors apporter un plus grand effet.

• 3. Si A est une matrice triangulaire supérieure ou inférieure, alors la méthode de substitution inverse est utilisée ou, par conséquent, la méthode de substitution directe, dans laquelle la composante de solution est trouvée à partir de la dernière (ou de la première équation), puis les composantes trouvées sont substitué dans les équations suivantes pour trouver les systèmes d'équations des composants de la solution suivante.
• 4. Si A - peut être réduit à une forme triangulaire par permutations, alors la réduction correspondante est effectuée, puis le système d'équations est résolu comme à l'étape 3.

5. Si A est une matrice symétrique (ou, dans le cas complexe, hermitienne) et que sa diagonale ne contient que des éléments réels positifs, alors selon que A est clairsemé ou non, l'élément a) ou l'élément b) est valable.

un. Si A n'est pas clairsemé, alors on tente d'effectuer la décomposition de Cholesky A = R T R suivie de la résolution de systèmes avec des matrices triangulaires R T et R : R T y = b et Rx = y. Dans ce cas, les procédures suivantes de la bibliothèque LAPACK sont appelées :

A et b sont réels et doubles - DLANGE, DPOTRF, DPOTRS, DPOCON
A et b sont complexes et doubles - ZLANGE, ZPOTRF, ZPOTRS, ZPOCON
A et b sont réels et de type unique - SLANGE, SPOTRF, SPOTRS, SDPOCON
Type complexe et unique A et B - CLANGE, CPOTRF, CPOTRS, CPOCON

La décomposition de Cholesky échouera si la matrice n'est pas définie positive. Mais il n'y a pas de contrôle préalable du caractère défini positif, et il est déterminé précisément au cours du processus de décomposition, car bien souvent il suffit d'effectuer plusieurs étapes de la décomposition de Cholesky pour clarifier le fait que la matrice n'est pas définie positive (dans le processus de calculs, il devient nécessaire d'extraire la racine carrée d'un nombre négatif). Si la décomposition de Cholesky n’est pas satisfaite, passez au point suivant.

b. Si A est clairsemé, alors les permutations symétriques des lignes et des colonnes sont d'abord effectuées à l'aide de l'algorithme de degré minimum symétrique (fonction symmmd) pour réduire le remplissage du facteur de décomposition de Cholesky, c'est-à-dire pour réduire le nombre de nouveaux éléments non nuls apparaissant lors du processus de remplissage :

p = symmmd(A) - le vecteur p contient des permutations

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

et deux systèmes avec multiplicateurs de Cholesky sont résolus, le premier avec un multiplicateur transposé et les permutations correspondantes dans le vecteur de droite

et le second, avec un facteur d'expansion et en entrant les composants de la solution dans les positions correspondantes dans le vecteur x
x(p, :) = R\y

• 6. Si A est une matrice de Hessenberg, alors elle est réduite à une matrice triangulaire supérieure (avec modification appropriée du membre de droite) puis le système est résolu par substitutions

7. Si A est une matrice carrée qui ne satisfait pas aux conditions des paragraphes 1 à 6, alors selon qu'elle est clairsemée ou non, les actions suivantes sont effectuées

un. Si A n'est pas une matrice clairsemée, alors par élimination gaussienne avec choix d'un élément leader (pour assurer la stabilité de la décomposition), une décomposition LU de la matrice A = LU est effectuée, où

U - matrice triangulaire supérieure
L est une matrice réduite à triangulaire par permutations de lignes

et la solution du système Ax = b est trouvée en résolvant séquentiellement des systèmes avec des matrices triangulaires Ly = b, Ux = y.
Pour effectuer la décomposition LU, les procédures suivantes de la bibliothèque LAPACK sont appelées :

A et b sont réels et doubles - DLANGE, DGESV, DGECON
A et b sont complexes et doubles - ZLANGE, ZGESV, ZGECON
A et b sont réels et de type unique - SLANGE, SGESV, SGECON
A et b sont de type complexe et unique - CLANGE, CGESV, CGECON

b. Si A est une matrice clairsemée, alors les colonnes sont réorganisées pour réduire le remplissage des facteurs L et U lors du processus de recherche de la décomposition. De plus, en réorganisant les lignes au cours du processus de décomposition LU, la stabilité de calcul est assurée, après quoi les systèmes avec des matrices triangulaires sont à nouveau résolus. Pour effectuer la décomposition LU d'une matrice clairsemée, les procédures de la bibliothèque UMFPACK sont utilisées

8. Si les points précédents de l'algorithme n'ont pas fonctionné et que, par conséquent, A n'est pas une matrice carrée, alors tout est déterminé par sa parcimonie et l'un des points ci-dessous fonctionne :

un. Si A n'est pas une matrice clairsemée, alors la décomposition QR AP = QR est effectuée, où P est la matrice de permutation de colonnes, Q est la matrice orthogonale (Q T Q = I) et R est triangulaire supérieure. Si A est de taille m sur n, alors Q est de taille m sur m et R est de taille m sur n. Ensuite, la solution du système se trouve comme suit :

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

puisque de Ax = b et AP = QR il résulte : QRx = bP et Rx = Q T bP.

b. Dans le cas d'une matrice A creuse et rectangulaire, cela ne sert à rien de calculer la matrice Q, puisqu'elle sera très probablement remplie. Par conséquent, l'algorithme de décomposition QR calcule c = Q T b (c'est-à-dire le membre de droite modifié) et résout le système Rx = c avec une matrice triangulaire pour obtenir une solution du système d'origine Ax = b.

En plus de tout ce qui précède, l'algorithme ci-dessus évalue la conditionnalité de la matrice et affiche un avertissement concernant une éventuelle erreur dans la solution si la matrice est mal conditionnée (voir section L'influence de la conditionnalité d'une matrice sur la précision de résoudre un système avec lui). L'algorithme ci-dessus contient diverses vérifications qui peuvent prendre un temps supplémentaire important. La section suivante, La fonction linsolve pour résoudre des systèmes d'équations linéaires, présente la fonction linsolve, qui vous permet de spécifier le type de matrice, réduisant ainsi le temps de calcul.

fonction linsolve pour résoudre des systèmes d'équations linéaires

Au lieu de résoudre un système (ou des systèmes avec plusieurs côtés droits définis dans une matrice) d'équations algébriques linéaires à l'aide du signe barre oblique inverse, vous pouvez utiliser la fonction linsolve, qui n'effectue pas toutes les vérifications matricielles inhérentes à l'algorithme d'opération \ ( voir section précédente).

La fonction linsolve appelée dans sa forme la plus simple avec deux arguments d'entrée et un argument de sortie
x = linsolve(A, b)
résout le système Ax = b de l'une des manières, selon que la matrice est carrée ou non.

• Si A est une matrice carrée, alors sa décomposition LU est d'abord calculée puis deux systèmes à matrices triangulaires L et U sont résolus.
• Si A est une matrice rectangulaire, alors sa décomposition QR est d'abord calculée puis le système à matrices triangulaires est résolu.

(Plus de détails sur les actions avec les facteurs d'expansion résultants sont écrits dans la section). De plus, si la matrice A est mal conditionnée ou si A est une matrice de rang incomplet, alors un avertissement correspondant s'affiche, par exemple :

UNE = ; b = ; x = linsolve(A,b) Avertissement : rang déficient, rang = 1, tol = 4.4686e-015. x = 0 0 2,0000

Pour supprimer ce message de l'impression dans la fenêtre de commande, appelez la fonction linsolve avec deux arguments de sortie
= linsolve(A, b)
alors la solution résultante s'écrira en x, et en r - soit le rang de la matrice (si A est une matrice rectangulaire), soit l'inverse de l'estimation de son numéro de condition (voir section), par exemple

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

Les principaux avantages de la fonction linsolve par rapport à l'opération \ apparaissent lors de la spécification du type de matrice du système d'équations. Pour ce faire, avant d'appeler la fonction linsolve, vous devez remplir une structure avec des informations sur la matrice avec les champs suivants

• SYM - symétrique ;
• LT - triangulaire inférieur ;
• UT - triangulaire supérieur ;
• UHESS-Hessenberg ;
• POSDEF - symétrique, défini positif ;
• RECT - rectangulaire ;
• TRANSA - s'il est nécessaire de résoudre un système avec une matrice transposée à une matrice donnée.

Chaque champ peut être vrai ou faux. Bien entendu, toutes les combinaisons de valeurs de champ ne sont pas valables ; par exemple, une matrice ne peut pas être triangulaire et en même temps symétrique et définie positive. Les combinaisons correctes de valeurs de champ sont rassemblées dans le tableau suivant

LT	Utah	UHESS	SYM	POSDEF	RECHERCHER	TRANSA
vrai	FAUX	FAUX	FAUX	FAUX	vrai faux	vrai faux
FAUX	vrai	FAUX	FAUX	FAUX	vrai faux	vrai faux
FAUX	FAUX	vrai	FAUX	FAUX	FAUX	vrai faux
FAUX	FAUX	FAUX	vrai	vrai	FAUX	vrai faux
FAUX	FAUX	FAUX	FAUX	FAUX	vrai faux	vrai faux

Si la matrice du système est définie positive, alors cela doit être pris en compte lors de la résolution, car pour les matrices définies positives, la solution est basée sur la décomposition de Cholesky, qui nécessite moins d'opérations que la décomposition LU utilisée lors de la résolution de systèmes avec des matrices carrées générales. . Ceci est facile à vérifier à l'aide de l'exemple suivant, dans lequel une matrice définie positive symétrique est créée (une matrice de nombres aléatoires est ajoutée avec une transposition et de grands nombres sont ajoutés à la diagonale) et un système de les équations avec cette matrice sont d'abord résolues comme un système avec une matrice de forme générale (opts. SYM = faux et opts.POSDEF = faux), puis comme avec une matrice définie symétrique et positive (opts.SYM = vrai et opts .POSDEF = vrai).

% définit tous les champs de la structure des opérations, sauf SYM et POSDEF opts.TRANSA = false ; opts.UHESS = faux; opts.RECT ​​​​= faux; opts.UT = faux ; opts.LT = faux ; % crée un vecteur de tailles de matrice N = 2.^(8:12); % crée des tableaux vides pour enregistrer les temps de solution TimeGeneral = ; TimePosSym = ; % dans une boucle nous créons des matrices et comparons les temps de solution pour n = N % nous créons une matrice définie symétrique et positive % et le vecteur du côté droit A = rand(n); A = A + A" + 2*n*eye(n); b = sum(A, 2); % résoudre le système comme un système avec une matrice générale opts.SYM = false; opts.POSDEF = false; Tstart = cputime; x = linsolve(A,b, opts); Tend = cputime; TimeGeneral = ; % résout le système comme un système avec une matrice de symétrie et de position opts.SYM = true; opts.POSDEF = true; Tstart = cputime; x = linsolve( A,b, opts); Tend = cputime; TimePosSym = ; end % affichage des graphiques temporels figure loglog(N, TimeGeneral, N, TimePosSym) legend("TimeGeneral", "TimePosSym")

Les coûts de calcul résultant de ces méthodes de résolution du système sont présentés dans le graphique ci-dessous

Bien entendu, si la matrice est triangulaire, alors il est très important de l'indiquer, puisque la résolution d'un système avec une matrice triangulaire s'effectue en opérations O(n 2), et si un algorithme de solution conçu pour une matrice de une forme générale est appliquée à un système à matrice triangulaire, il faudra alors environ O(n 3) opérations.

La fonction linsolve ne vérifie pas si la matrice satisfait aux propriétés spécifiées, ce qui peut entraîner une erreur. Par exemple, si, lors de la résolution d'un système avec la matrice suivante et le membre de droite

UNE = ; b = ; définissez true pour le champ UT (et définissez tous les autres sur false) opts.UT = true; opts.TRANSA = faux; opts.LT = faux ; opts.UHESS = faux; opts.SYM = faux ; opts.POSDEF = faux ; opts.RECT ​​​​= faux; puis lors de la résolution du système, la fonction linsolve considérera la matrice comme une matrice triangulaire supérieure et sélectionnera les éléments dont elle a besoin dans le triangle supérieur A x = linsolve(A,b, opts) Cela donnera une solution x = 1 1 1 correspondant à la matrice A = ;

L'influence de la conditionnalité matricielle sur la précision de la résolution d'un système avec celle-ci

Dans cette section, nous examinons comment les erreurs dans les éléments du vecteur de droite et dans la matrice d’un système d’équations linéaires Ax = b peuvent affecter la solution de ce système. Passons d'abord au cas de l'introduction de perturbations dans le vecteur du membre droit b. Nous résolvons donc deux systèmes

De plus, on suppose que nous résolvons exactement chacun des systèmes, et la principale question qui se pose est de savoir dans quelle mesure la solution x du système (1) sera différente de la solution du système (2) avec une perturbation sur le côté droit δb. C'est une question assez importante, puisque les éléments du côté droit peuvent être obtenus à la suite de certaines mesures, contenant respectivement l'erreur δb k dans chaque composante b k. J'aimerais qu'en mesurant la même quantité (à chaque fois avec ses propres petites erreurs), les solutions correspondantes de systèmes avec des côtés droits légèrement différents ne diffèrent pas non plus beaucoup les unes des autres. Malheureusement, ce n'est pas toujours le cas. La raison en est une caractéristique de la matrice appelée son conditionnalité. Ceci sera discuté plus loin.

Tout d'abord, vous devez introduire une mesure de la proximité des vecteurs et de x, c'est-à-dire vecteur d'erreur. La mesure de la grandeur d'un vecteur est la norme (elle peut être définie de différentes manières). Pour l'instant, prenons comme norme la norme euclidienne habituelle d'un vecteur (la racine carrée de la somme des carrés de ses composantes), c'est-à-dire

Par conséquent

où n est le nombre d'inconnues et d'équations dans le système. En plus de la norme vectorielle pour estimer la norme du vecteur, nous avons également besoin d’une norme matricielle pour estimer la norme de la matrice. Prenons la norme matricielle, qui est définie comme suit (on l’appelle la norme spectrale) :

ceux. la norme spectrale de la matrice est égale à la racine carrée de la valeur propre maximale de la matrice A T A. MATLAB dispose d'une fonction de norme pour calculer les normes des matrices et des vecteurs, qui permet notamment de calculer les normes ci-dessus. La raison pour laquelle nous avons choisi ces normes particulières de vecteurs et de matrices est expliquée en détail dans la section, où quelques calculs et définitions sont donnés. Ceci est lié à l'estimation que nous utiliserons pour l'erreur dans la résolution du système (la dérivation de cette inégalité est également donnée dans la section mentionnée) :

Nous désignons ici le numéro de condition de la matrice, qui est défini comme suit :

De l'inégalité ci-dessus, il s'ensuit que plus le numéro de condition de la matrice système est grand, plus l'erreur relative dans la solution peut être grande avec une petite erreur du côté droit.

Considérons un exemple classique de matrice mal conditionnée - la matrice de Hilbert - et découvrons comment l'erreur du côté droit du système affecte l'erreur dans la solution. La matrice de Hilbert est définie comme suit

Pour créer une matrice de Hilbert, MATLAB fournit la fonction hilb, dont l'argument d'entrée spécifie la taille de la matrice. Prenons une petite matrice 6 par 6 :

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,1 111 0,1000 0,0909

X = uns(n, 1); b = H*x b = 2,4500 1,5929 1,2179 0,9956 0,8456 0,7365

On voit qu'il n'y a rien de « suspect » ni dans la matrice ni dans le vecteur de droite ; tous les nombres ne sont pas très différents les uns des autres. Formons maintenant le membre droit perturbé b + δb en ajoutant des petits nombres de l'ordre de 10 -5 au vecteur b et résolvons le système avec le membre droit perturbé pour obtenir le vecteur .

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

On peut voir que la solution résultante est loin d’être celle exacte où toutes les unités devraient se trouver. Calculons l'erreur relative dans la solution et sur le côté droit (la fonction norme calcule par défaut la norme euclidienne du vecteur) :

Delta_x = x - x_tilda ; GAUCHE = norme(delta_x)/norme(x) GAUCHE = 1,1475 DROITE = norme(delta_b)/norme(b) DROITE = 2,7231e-005

Ainsi, l’erreur dans la solution est de l’ordre de l’unité, bien que les changements du côté droit soient de l’ordre de 10 -5. Cela correspond parfaitement à l’inégalité d’erreur ci-dessus. En effet, calculons le numéro de condition cond(H) à l'aide de la fonction MATLAB appelée cond et calcule par défaut pour la norme spectrale de la matrice

C = cond(H) c = 1,4951e+007

GAUCHE ≈ 1.1475 moins

C*N'est-ce pas ? 1,4951e+07 * 2,7231e-05 ≈ 407

et l'inégalité est satisfaite (même avec une certaine marge).

À mesure que la taille de la matrice de Hilbert augmente, l'erreur dans la solution ne fera qu'augmenter (cela est facile à vérifier en définissant n = 7, 8, ...). De plus, pour n = 12 un message s'affichera indiquant que la matrice est mal conditionnée et que la solution est peut-être incorrecte

Attention : La matrice est proche du singulier ou mal mise à l'échelle.
Les résultats peuvent être inexacts. RCOND = 2.409320e-017.

Comme mesure de conditionnalité, la valeur choisie ici est RCOND égale à un divisé par l'estimation du numéro de condition (le numéro de condition est estimé à l'aide d'un algorithme assez rapide qui fonctionne beaucoup plus vite que cond, qui est décrit plus en détail dans l'aide sur la fonction recnd). Si la valeur de RCOND est petite, alors la matrice est considérée comme mal conditionnée.

L'augmentation de l'erreur dans la solution à mesure que la taille de la matrice de Hilbert augmente est due au fait que le numéro de condition de la matrice de Hilbert augmente très rapidement avec sa taille. Ceci est facile à vérifier à l'aide d'une simple boucle et de la fonction sémilogie (l'échelle le long de l'axe des ordonnées est logarithmique) :

N = 1:20 ; C = zéros (1, 20) ; pour n = N H = hilb(n); C(n) = cond(H); fin de la grille de sémilogie(N, C) activée, titre("cond(H)"), xlabel("n")

À propos, puisque la fonction cond trouve le numéro de condition à l'aide d'une méthode numérique (à savoir, le développement singulier pour trouver des nombres singuliers), le numéro de condition après n = 12 n'est plus calculé correctement ; en fait, il devrait encore croître, ce qui peut être vérifié à l'aide de calculs symboliques dans MATLAB et d'opérations avec un nombre donné de chiffres significatifs

N = 1:20 ; C = zéros(1, N); chiffres (60) pour n = N H = vpa(sym(hilb(n))); % de calcul de la matrice de Hilbert au 60ème chiffre près sigma = svd(H) ; % trouver les valeurs singulières de la matrice de Hilbert % calculer le numéro de condition de la matrice de Hilbert C(n) = max(double(sigma))/min(double(sigma)); fin de la grille de sémilogie(N, C) activée, titre("cond(H)"), xlabel("n")

Considérons maintenant trois points importants concernant cet exemple.

Le premier d'entre eux - résoudre le système Hx = b (avec le vecteur de droite correspondant à la solution à partir des unités) à l'aide de l'opérateur barre oblique inverse ne donnera pas les unités exactes, l'erreur sera déjà dans le dixième chiffre (bien que dans MATLAB tout les calculs sont effectués avec une double précision par défaut)

Format long H\b ans = 0,999999999999916 1,00000000002391 0,99999999983793 1,00000000042209 0,999999999953366 1,00000000018389

Cela se produit parce que lors du calcul du vecteur b en multipliant la matrice H par un vecteur de tous les uns, nous y avons déjà inclus une erreur. De plus, les erreurs d'arrondi dans le processus de résolution du système ont également joué un rôle et une mauvaise conditionnalité de la matrice (même de petite taille) a conduit à de telles erreurs dans la solution. En effet, pour les petites matrices avec un petit numéro de condition, un tel effet ne sera pas observé. Prenons une matrice 6 x 6 de nombres aléatoires et calculons son numéro de condition

A = rand(n); cond(A) ans = 57,35245279907571

Ensuite on forme le membre de droite correspondant à la solution exacte de toutes les unités

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

et résoudre le système, ce qui donne une bonne précision

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

Le deuxième point important est lié au déterminant de la matrice. Calculons le déterminant de la matrice de Hilbert 6 x 6 avec laquelle nous avons résolu le système

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

Le déterminant s'avère suffisamment petit, à partir duquel on peut conclure à tort que la source du problème d'une grande erreur dans la solution avec une petite perturbation du côté droit du système est la petitesse du déterminant. En réalité, ce n'est pas le cas ; si le déterminant du système est petit, cela ne veut pas dire que la matrice est mal conditionnée et, par conséquent, il y a une grande erreur dans la solution. En effet, prenons la matrice

UNE = ;

Son déterminant est très petit, environ 10 -121

Dét(A) ans = 9,9970e-121

(le mot « très » est bien sûr conditionnel, mais il est plus petit que le déterminant d'une matrice de Hilbert 6 x 6, avec laquelle nous avons eu des problèmes lors de la résolution du système). Cependant, la petitesse du déterminant n'affectera en rien l'erreur dans la solution lorsque le côté droit du système est perturbé, ce qui est facile à montrer en formant un système avec une solution connue, en introduisant des perturbations dans le côté droit et en résolvant le système:

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

Ainsi, l'erreur relative dans la solution est presque égale à l'erreur relative du côté droit, puisque le numéro de condition de la matrice A ci-dessus est légèrement supérieur à 1, à savoir :

C = cond(A) c = 1,0303

Et enfin, considérons la troisième question concernant l'obtention du signe égal dans l'inégalité pour une erreur dans la solution

Autrement dit, découvrons : peut-il y avoir une situation où nous produisons une petite perturbation relative du côté droit du système, disons 10 -5, le numéro de condition de la matrice du système est égal à 10 10 , et dans la solution, une erreur relative de 10 5 est obtenue. Vous pouvez parcourir de nombreux exemples, essayer différentes matrices, ne pas atteindre l'égalité et conclure à tort qu'il s'agit simplement d'une surestimation d'en haut pour une erreur dans la solution. Cependant, ce n’est pas le cas, comme l’exemple suivant nous en convainc.

dans lequel la perturbation relative du côté droit est de 10 -5

DROITE = norme(delta_b)/norme(b) DROITE = 1.0000e-005

L'erreur relative dans la résolution du système s'avère être 10 5

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

Et cela arrive parce que
1) dans cet exemple, le numéro de condition de la matrice A est 10 10 ;

C = cond(A) c = 1,0000e+010

2) l'inégalité pour l'erreur dans la solution se transforme en égalité.

Si nous définissons le format sur long, nous voyons que LEFT, RIGHT et le numéro de condition ne sont pas exactement 10 5 , 10 -5 et 10 10 , respectivement, mais cela est dû à des erreurs d'arrondi. Si elle était résolue en arithmétique exacte, l'égalité serait exacte, ce qui peut être montré analytiquement (voir section).

Dans les exemples précédents, nous avons supposé que la matrice ne contenait pas d'erreurs et que les erreurs ne pouvaient se trouver que sur le côté droit du système d'équations linéaires. En pratique, il peut s'avérer que les éléments de la matrice système sont spécifiés avec des erreurs, c'est-à-dire au lieu du système

|| Hache || V = || b || V ⇒ || Un || M || x || V ≥ || b || V

où || || M est une norme matricielle. Pour ce faire, la norme matricielle || || M et norme vectorielle || || V doit satisfaire l'inégalité suivante

|| Un || M || x || V ≥ || Hache || V

pour toutes matrices A et vecteurs x. Si cette inégalité est vraie, alors la norme matricielle || || M s'appelle convenu avec norme vectorielle || || V. On sait par exemple que la norme spectrale

(racine carrée de la valeur propre maximale de la matrice A T) est conforme à la norme vectorielle euclidienne

C'est pourquoi dans cette section nous avons réalisé toutes les expériences en utilisant ces normes.

Ayant reçu l'inégalité || Un || M || x || V ≥ || b || V notons en outre que de Ax = b il s'ensuit que . Puisque la matrice A est non singulière, il s’ensuit que δx = A -1 δb et || δx || V = || A -1 δb || V. Encore une fois, nous utilisons la cohérence des normes et obtenons l'inégalité || δx || V ≤ || A-1 || M || δb || V. De plus, dans les deux inégalités obtenues

|| Un || M || x || V ≥ || b || V et || δx || V ≤ || A-1 || M || δb || V

Nous divisons la plus petite valeur de l'une des inégalités par la plus grande valeur de l'autre inégalité, et la plus grande valeur, respectivement, par la plus petite.

et par une simple transformation on obtient finalement l'inégalité recherchée

où cond(A) = || Un || M* || A-1 || M.

Notez que lors du calcul de l’inégalité pour l’erreur, nous avons utilisé uniquement le fait que la norme matricielle est cohérente avec la norme vectorielle. Cela est vrai non seulement pour la norme spectrale d’une matrice et la norme euclidienne d’un vecteur, mais aussi pour d’autres normes. Ainsi, par exemple, la norme maximale de la matrice de colonnes, calculée par la formule

cohérent avec la première norme vectorielle

et la norme maximale de la matrice de lignes, calculée par la formule

cohérent avec la norme vectorielle infinie

La fonction norme calcule non seulement la norme euclidienne d'un vecteur et la norme spectrale d'une matrice, mais également les normes vectorielles et matricielles répertoriées ci-dessus. Pour ce faire, vous devez l'appeler avec un deuxième argument supplémentaire :

• q = norm(A, 1) - norme de colonne maximale de la matrice A
• q = norm(A, inf) - norme de ligne maximale de la matrice A
• q = norm(a, 1) - première norme vectorielle de a
• q = norm(a, inf) - norme vectorielle infinie a

Numéro de condition matricielle cond(A) = || Un || M* || A-1 || M par rapport à diverses normes matricielles peut être calculé à l'aide de la fonction cond. Si cond est appelé avec un argument d'entrée (une matrice), alors le numéro de condition est calculé par rapport à la norme de la matrice spectrale. L'appel de la fonction cond avec un argument supplémentaire renvoie les numéros de condition relatifs à la norme matricielle spécifiée :

• c = cond(A, 1) - numéro de condition par rapport à la norme de colonne maximale de la matrice
• с = cond(A, inf) - numéro de condition par rapport à la norme de ligne maximale de la matrice

A titre d'exemple démontrant l'importance de la cohérence de la norme matricielle avec la norme vectorielle dans l'inégalité d'erreur, nous donnons un exemple avec la matrice A, le vecteur de droite b et le vecteur d'erreur du côté droit delta_b.

UNE = ; 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];

Calculons les parties droite et gauche de l'estimation en utilisant la première norme vectorielle, et le numéro de condition de la matrice par rapport à la norme de la matrice spectrale

DROITE = norme(delta_b,1)/norme(b,1); c = cond(UNE); x = A\b; x_tilda = A\(b+delta_b); delta_x = x - x_tilda ; GAUCHE = norme(delta_x,1)/norme(x,1); formater court et disp()

On obtient les valeurs suivantes pour les côtés gauche et droit de l'inégalité
9.9484e+004 9.9323e+004

Le côté gauche ne doit pas dépasser le droit (comme prouvé ci-dessus), mais il s'est avéré plus grand en raison du choix de normes incohérentes.

Examinons maintenant pourquoi le numéro de condition d'une matrice ne peut pas être inférieur à un. Le numéro de condition de la matrice A est défini comme cond(A) = || Un || M* || A-1 || M , où || Un || M est une norme matricielle de A. La norme matricielle (c'est-à-dire la règle par laquelle un nombre est associé à chaque matrice) ne peut pas être arbitraire ; elle doit satisfaire les quatre axiomes suivants :

A1. || Un || ≥ 0 pour toute matrice A et || Un || = 0 si et seulement si A est une matrice nulle.

A2. || αA || = | α | * || Un || pour toute matrice A et nombre α.

A3. || A + B || ≤ || Un || + || B || pour toutes les matrices A et B

A4. || AB || ≤ || Un || * || B || pour toutes les matrices A et B.

D'après le dernier axiome, il est clair que la norme n'est définie que pour les matrices carrées (bien que dans les formules ci-dessus pour calculer diverses normes, en principe, une telle restriction n'existe pas). De plus, du dernier axiome il résulte que toute norme de la matrice identité I n'est pas inférieure à un, voire

|| Je || = || Je*je || ≤ || Je || 2 ⇒ || Je || ≥ 1.

Puis, toujours en utilisant le quatrième axiome, on constate que le numéro de condition de la matrice est toujours supérieur à un (vrai pour le numéro de condition de la matrice par rapport à une norme matricielle arbitraire)

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

Complétons cette section en examinant la raison de l'apparition d'une égalité exacte dans l'estimation de l'erreur de la solution lorsque le côté droit est perturbé :

et créer des exemples correspondants dans MATLAB. (La discussion ci-dessous est contenue, par exemple, dans le livre de J. Forsythe, K. Mohler. Solution numérique des systèmes d'équations algébriques linéaires. M : « Mir », 1969.)

Un rôle essentiel dans le raisonnement est joué par le théorème de décomposition singulière d'une matrice, selon lequel pour toute matrice réelle de taille n il existe deux matrices orthogonales U et V de taille n par n (U T U=UU T et V T V = VV T) tel que le produit D = U T AV soit une matrice diagonale, et on peut choisir U et V tels que

où μ 1 ≥ μ 2 ≥…≥μ r ≥ μ r+1 =…=μ n =0,

et r est le rang de la matrice A. Les nombres μ k sont appelés les nombres spectraux de la matrice A. Pour une matrice non singulière A, ce qui suit est vrai :

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

Le prochain fait important est que la multiplication par une matrice orthogonale ne change pas la norme euclidienne du vecteur, c'est-à-dire pour tout vecteur x avec n éléments et toute matrice orthogonale U de taille n par n l'égalité est vraie

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

Puisque la multiplication par une matrice orthogonale ne change pas la norme spectrale, alors

par conséquent, pour le numéro de condition de la matrice, l’égalité suivante est vraie :

Nous avons deux systèmes : Ax = b (avec solution exacte x) et (avec solution exacte ). Évidemment, l’erreur δx est une solution à un système dont le côté droit est une perturbation δb, c’est-à-dire systèmes Aδx = δb. Soit D = U T AV la décomposition en valeurs singulières de la matrice A, alors UDV T = A du fait que U et V sont des matrices orthogonales. Plus loin

Ax = b ⇒ UDV T x = b.

Introduisons la notation

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

alors les systèmes suivants sont équivalents

Ax = b ⇔ Dx" = b"

De manière tout à fait similaire, considérons le système par rapport à l'erreur

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

Nous introduisons la notation

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

pour lequel les systèmes suivants s'avèrent équivalents

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

Ceux. nous avons obtenu des systèmes équivalents simples avec une matrice diagonale, sur la diagonale de laquelle se trouvent les nombres spectraux de la matrice A.

Choisissons maintenant les membres droits de ces systèmes d'une manière particulière, à savoir :

où β > 0, alors la solution correspondante au système Dx" = b" sera

où μ 1 est la valeur singulière maximale de la matrice A. On fera de même avec le système Dδx" = δ b", à savoir soit

où γ > 0, alors la solution correspondante au système Dδx" = δb" sera

De la préservation de la norme d'un vecteur lorsqu'il est multiplié par une matrice orthogonale, il s'ensuit que

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

Exactement de la même manière on obtient les égalités :

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

et depuis

puis nous obtenons enfin

Ainsi, pour chaque matrice A, il est possible de construire un vecteur sur le côté droit du système et sa perturbation tel que l'erreur dans la solution sera le produit du numéro de condition de la matrice et de l'erreur sur le côté droit. Dans MATLAB, la construction correspondante peut être effectuée sans utiliser de décomposition en valeurs singulières (bien qu'elle puisse être trouvée en utilisant la fonction svd).

Nous définissons d’abord n ​​et obtenons deux matrices orthogonales U et V en effectuant une décomposition QR de matrices n par n à partir de nombres aléatoires :

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

D = diag();

Après cela, nous construisons la matrice A en utilisant la matrice diagonale D et les matrices orthogonales U et V

A = U*D*V" ;

Le numéro de condition de la matrice A coïncidera avec le numéro de condition de la matrice D, qui dans notre exemple est égal à 10 10

Bêta = 1 ; gamma = 1e-5 ;

et construire les vecteurs b" et δb"

B1 = "; db1 = ";

à partir duquel on retrouve les vecteurs b et δb

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

calculer les côtés gauche et droit de l'inégalité

dx = x - x_tilda ; DROITE = norme(db)/norme(b); GAUCHE = norme(dx)/norme(x);

et fais-les sortir

Formater court et disp()

Nous obtenons l'égalité

1.Fenêtre de commande(Fenêtre de commande).

Les expressions mathématiques sont écrites sur la ligne de commande après l'invite " >> ". Par exemple,

Pour effectuer l'action, appuyez sur la touche « Entrée ».

Par défaut, le programme écrit le résultat dans la variable spéciale ans.

Pour enregistrer le résultat sous un autre nom, utilisez le nom de la variable et le signe d'affectation « = », par exemple

>> z = 1,25 /3,11

Modifier dans Fenêtre de commande Vous ne pouvez utiliser que la ligne de commande actuelle. Pour éditer une commande saisie précédemment, vous devez placer le curseur dans la ligne de saisie et utiliser les touches « » ou « ».

Si une commande se termine par « ; », alors le résultat de son action n’est pas affiché sur la ligne de commande.

La fenêtre de commande peut être fermée avec le bouton « » et le bouton « » sert à séparer la fenêtre de l'interface système. Vous pouvez revenir au formulaire de fenêtre standard en utilisant le menu :

Menu principal → Bureau → Disposition du bureau → Défaut.

Vous pouvez effacer la fenêtre de commande en utilisant le menu :

Menu principal → Modifier → Effacer la fenêtre de commande.

Menu principal du système MatLab.

Menu principal MatLab contient les éléments suivants :

Déposer(Fichier) – travailler avec des fichiers ;

Modifier(Modifier) ​​– édition ;

Voir(Affichage) – gestion des fenêtres ;

la toile– communication avec la société de développement via Internet ;

Fenêtre(Fenêtre) – connexion avec les fenêtres du système ;

Aide(Aide) – lien vers le système d'aide MatLab.

Barre d'outils du système MATLAB.

Les boutons de la barre d'outils ont les objectifs suivants :

Nouveau fichier(Nouveau) – affiche les fenêtres de l'éditeur de fichiers ;

Fichier ouvert(Ouvrir) – ouvre les fenêtres de téléchargement de fichiers ;

Couper(Couper) – coupe le fragment sélectionné et le place dans le presse-papiers ;

Copie(Copier) – copie le fragment sélectionné dans le presse-papiers ;

Pâte(Coller) – transfère le fragment sélectionné du presse-papiers vers la ligne d'entrée ;

annuler(Annuler) – annule le résultat de l'opération précédente ;

Refaire(Répéter) – répète le résultat de l'opération précédente ;

Simulink– crée un modèle Simulink (extensions MatLab);

Fenêtre d'aide(Aide) – ouvre les fenêtres d'aide.

4. Format de sortie des résultats de calcul .

Lors de la saisie de nombres réels, un point est utilisé pour séparer la partie fractionnaire!

>> s = 0,3467819

Le résultat du calcul est affiché au format court(notation courte pour un nombre), qui est définie par défaut. Vous pouvez modifier le format en long(notation de nombres longs).

>> format long

0.34678190000000

Sur la liste Format numérique formats de nombres disponibles

court– une brève mention du numéro ;

long– notation des nombres longs ;

court e– notation courte d'un nombre au format virgule flottante ;

long e– enregistrement long d'un nombre au format virgule flottante ;

court g– la deuxième forme d'une notation courte d'un nombre ;

longue g– la deuxième forme de notation longue d'un nombre ;

Le format d'affichage des données numériques peut être défini dans le menu Déposer(fichier) élément Préférences(préférences). Aller à l'onglet Fenêtre de commande(fenêtre de commande). En option Affichage du texte(affichage de texte), dans la liste Format numérique(format numérique) défini court g, en option Affichage numérique(afficher les numéros) compact. Ces formats de sortie génèrent des nombres sous une forme universelle de cinq chiffres significatifs et suppriment les espaces blancs entre les lignes.

Bases des calculs dans MatLab.

Pour effectuer des opérations arithmétiques simples dans MatLab les opérateurs sont utilisés :

· addition et soustraction +, – ;

· Multiplication et division *, / ;

· exponentiation ^ .

Quelques variables spéciales:

ans – le résultat de la dernière opération sans signe d'affectation ;

eps – erreur relative dans les calculs en virgule flottante ;

pi – nombre ;

je ou j – unité imaginaire ;

Inf – l'infini ;

NaN – valeur indéfinie.

Quelques fonctions élémentaires intégréesMatLab:

exp(x) – exposant du nombre x ;

log(x) – logarithme népérien du nombre x ;

sqrt(x) – racine carrée du nombre x ;

abs(x) – module du nombre x ;

sin(x), cos(x), tan(x), cot(x) – sinus, cosinus, tangente, cotangente du nombre x ;

asin(x), acos(x), atan(x), acot(x) – arc sinus, arccosinus, arctangente, arccotangente du nombre x ;

sec(x), csc(x) – sécante, cosécante du nombre x ;

round(x) – arrondir le nombre x à l'entier le plus proche ;

mod(x,y) – reste de la division entière de x par y, en tenant compte du signe ;

sign(x) – renvoie le signe du nombre x.

Calculons la valeur de l'expression

>> exp(–2,5)*log(11,3)^0,3 – sqrt((sin(2,45*pi)+cos(3,78*pi))/tan(3,3))

Si l'opérateur ne peut pas être placé sur une ligne, alors il est possible de continuer à le saisir sur la ligne suivante si vous indiquez le signe de suite « … » à la fin de la première ligne, par exemple,

>> exp(–2,5)*log(11,3)^0,3 – ...

carré((sin(2,45*pi)+cos(3,78*pi))/tan(3,3))

Fonctions pour travailler avec des nombres complexes:

Saisir un nombre complexe

>> z = 3 + 4i

3.0000 + 4.0000i

Les fonctions abs(z), angle(z) renvoient le module et l'argument d'un nombre complexe, où , ;

complex(a,b) construit un nombre complexe à partir de ses parties réelles et imaginaires :

conj(z)renvoie le nombre conjugué complexe ;

imag(z), real(z) renvoie la partie imaginaire et réelle du nombre complexe z.

6. Vecteurs.

Saisir, additionner, soustraire, multiplier par un nombre.

Vecteur dans MatLab formé à l’aide de l’opérateur crochets. Dans ce cas, les éléments d'un vecteur colonne sont séparés par un point-virgule « ; », et les éléments d'un vecteur ligne sont séparés par un espace « » ou une virgule « »,.

Introduisons un vecteur colonne.

>> x =

Introduisons un vecteur ligne .

>> y =

Pour transposer un vecteur, utilisez l'apostrophe « ' » :

>> z = y'

Pour trouver la somme et la différence des vecteurs, utilisez les signes « + » et « – » :

>> c = x + z

La multiplication d'un vecteur par un nombre s'effectue aussi bien à droite qu'à gauche à l'aide du signe « * ».

Les vecteurs peuvent être des arguments pour des fonctions intégrées, par ex.

>> d = péché(c)

Pour désigner les éléments des vecteurs, des parenthèses () sont utilisées, par exemple,

>> x_2 = x(2)

Le dernier élément du vecteur peut être sélectionné en tapant la commande

>> X_end = x(fin)

A partir de plusieurs vecteurs vous pouvez en créer un, par exemple

>> r =

1.3 5.4 6.9 7.1 3.5 8.2

Le caractère deux-points " : " est utilisé pour séparer plusieurs éléments d'un vecteur, par exemple

>> w = r(3:5)

Le caractère deux-points " : " permet également de remplacer des éléments d'un vecteur, par exemple,

>> r(3:5)= 0

1.3 5.4 0 0 0 8.2

Le symbole « : » peut également être utilisé pour construire un vecteur dont chaque élément diffère du précédent par un nombre constant, c'est-à-dire étape, par exemple

>> h =

1 1.2 1.4 1.6 1.8 2

Le pas peut être négatif (dans ce cas, le nombre de départ doit être supérieur au nombre final).

Un pas égal à un peut être omis

>> k =

Fonctions de base pour travailler avec des vecteurs.

• length(x) – déterminer la longueur du vecteur x ;
• prod(x) – multiplication de tous les éléments du vecteur x ;
• sum(x) – somme de tous les éléments du vecteur x ;
• max(x) – trouver l'élément maximum du vecteur x ;
• min(x) – trouver l’élément minimum du vecteur x.

Si vous appelez la fonction min ou max avec deux arguments de sortie = min(x),

puis la première variable se voit attribuer la valeur de l'élément minimum (maximum), et la deuxième variable se voit attribuer le numéro de cet élément.

7 matrices.

Diverses méthodes de saisie matricielle.

1. Une matrice peut être saisie sous forme de vecteur colonne composé de deux éléments, chacun étant un vecteur ligne et séparé par un point-virgule. Par exemple, introduisons la matrice

>> A =

2. La matrice peut être saisie ligne par ligne en exécutant la séquence de commandes :

>> A = )