Des articles

4 : Comportement à proximité de trajectoires - Linéarisation - Mathématiques


Nous allons maintenant discuter d'une méthode d'analyse de la stabilité qui utilise la linéarisation sur l'objet dont la stabilité est d'intérêt. Pour l'instant, les « objets d'intérêt » sont des solutions spécifiques d'un champ vectoriel. La structure des solutions de systèmes linéaires à coefficients constants est traitée dans de nombreux manuels ODE. Le livre d'Arnold est également très bon, mais la présentation est plus compacte, avec moins d'exemples.

On commence par considérer un champ de vecteurs général non autonome :

[dot{x} = f(x, t), x in mathbb{R}^n, label{4.1}]

et on suppose que

[ar{x}(t, t_{0}, x_{0}), label{4.2}]

est la solution de l'équation ef{4.1} pour laquelle on souhaite déterminer ses propriétés de stabilité. Comme lorsque nous avons introduit les définitions de stabilité, nous procédons en localisant le champ de vecteurs autour de la solution d'intérêt. Nous le faisons en introduisant le changement de coordonnées

(x = y+ar{x})

pour l'équation ef{4.1} comme suit :

(dot{x} = dot{y}+dot{ar{x}} = f(y+ar{x}, t)),

ou alors

(dot{y} = f(y+ar{x}, t)dot{ar{x}}),

[= f(y+ar{x}, t)f(ar{x}, t), label{4.3}]

où nous omettons les arguments de (ar{x}(t, t_{0}, x_{0})) pour une notation moins lourde. Ensuite, nous développons (f(y+ar{x}, t)) en (y) à propos de la solution (ar{x}), mais nous n'exigerons explicitement que les termes d'ordre principal

[f(y+ar{x}, t) = f(ar{x}, t)+Df(ar{x},t)y+mathbb{O}(|y|^2), étiquette{4.4}]

où (Df) désigne la dérivée (c'est-à-dire la matrice Jacobienne) de la fonction à valeur vectorielle f et (mathbb{O}(|y|^2)) désigne des termes d'ordre supérieur dans le développement de Taylor dont nous n'aurons pas besoin sous forme explicite. En substituant ceci dans l'équation ef{4.4} donne :

(dot{y} = f(y+ar{x}, t) - f(ar{x}, t)),

(= f(ar{x}, t)+Df(ar{x}, t)y+mathbb{O}(|y|^2)f(ar{x}, t)),

[= Df(ar{x}, t)y+mathbb{O}(|y|^2). label{4.5}]

Gardez à l'esprit que nous nous intéressons au comportement des solutions proches de (ar{x}(t, t_{0}, x_{0})), c'est-à-dire pour (y) petit. Par conséquent, dans cette situation, il semble raisonnable que négliger le (mathbb{O}(|y|^2)) dans l'équation ef{4.5} serait une approximation qui nous fournirait l'information particulière que nous recherchons. Par exemple, nous fournirait-il suffisamment d'informations pour déterminer la stabilité ? En particulier,

[dot{y} = Df(ar{x}, t)y, label{4.6}]

est appelée la linéarisation du champ vectoriel (dot{x} = f (x, t)) autour de la solution (ar{x}(t, t_{0}, x_{0}) ).

Avant de répondre à la question de savoir si l'équation ef{4.1} fournit ou non une approximation adéquate des solutions de l'équation ef{4.5} pour y "petit", nous allons d'abord étudier les champs de vecteurs linéaires seuls.

Les champs de vecteurs linéaires peuvent également être classés comme non autonome ou alors autonome. Les champs de vecteurs linéaires non autonomes sont obtenus en linéarisant un champ de vecteurs non autonomes autour d'une solution (et en ne retenant que les termes linéaires). Ils ont la forme générale :

[dot{y} = A(t)y, y(0) = y_{0}, label{4.7}]

[A(t) equiv Df(ar{x}(t, t_{0}, x_{0}), t) label{4.8}]

est une matrice (n imes n). Ils peuvent également être obtenus en linéarisant un champ de vecteurs autonome autour d'une solution dépendante du temps.

Un champ vectoriel linéaire autonome est obtenu en linéarisant un champ vectoriel autonome autour d'un point d'équilibre. Plus précisément, soit (dot{x} = f(x)) un champ de vecteurs autonome et soit (x = x_{0}) un point d'équilibre, soit (f(x_{0}) = 0). Le champ vectoriel autonome linéarisé autour de ce point d'équilibre a la forme :

[dot{y} = Df(x_{0})y, y(0) = y_{0}, label{4.9}]

ou alors

[dot{y} = Ay, y(0) = y_{0}, label{4.10}]

où (A equiv Df(x_{0})) est une matrice (n imes n) de nombres réels. Ceci est important car (4.10) peut être résolu en utilisant des techniques d'algèbre linéaire, mais l'équation ef{4.7}, en général, ne peut pas être résolue de cette manière. Par conséquent, nous allons maintenant décrire la solution générale de (4.10).

La solution générale de l'équation ef{4.10} est donnée par :

[y(t) = e^{At}y_{0}. label{4.11}]

Afin de vérifier que c'est la solution, il suffit de substituer le côté droit et le côté gauche de (4.10) et de montrer que l'égalité est vérifiée. Cependant, nous devons d'abord expliquer ce qu'est (e^{At}), c'est-à-dire l'exponentielle de la matrice (n imes n) A (en examinant l'équation ef{4.11} il devrait être clair que si l'équation ef{4.11} doit avoir un sens mathématique, alors (e^{At}) doit être une matrice (n fois n)).

Tout comme l'exponentielle d'un scalaire, l'exponentielle d'une matrice est définie par la série exponentielle comme suit :

(e^{At} equiv mathbb{I}+At+frac{1}{2!}A^{2}t^{2}+···+frac{1}{n!}A ^{n}t^{n}+cdots),

[= sum_{i=0}^{n}frac{1}{i!}A^{i}t^{i}, label{4.12}]

où (mathbb{I}) désigne la matrice identité (n imes n). Mais nous devons encore répondre à la question : « cette série exponentielle impliquant des produits de matrices a-t-elle un sens mathématique » ? Certes, nous pouvons calculer des produits de matrices et les multiplier par des scalaires. Mais nous devons donner un sens à une somme infinie de tels objets mathématiques. Nous faisons cela en définissant la norme d'une matrice puis en considérant la convergence de la série en norme. Lorsque cela est fait, le "problème de convergence" est exactement le même que celui de l'exponentielle d'un scalaire. Par conséquent, la série exponentielle pour une matrice converge absolument pour tout t, et donc elle peut être différenciée par rapport à t terme par terme, et la série résultante de dérivées converge également de manière absolue.

Ensuite, nous devons soutenir que l'équation ef{4.11} est une solution de l'équation ef{4.10}. Si on différencie la série (4.12) terme à terme, on obtient que :

[frac{d}{dt}e^{At} = Ae^{At} = e^{At}A, label{4.13}]

où nous avons utilisé le fait que les matrices A et (e^{At}) commutent (cela se déduit facilement du fait que A commute avec n'importe quelle puissance de A) Il résulte alors de ce calcul que :

[dot{y} = frac{d}{dt}e^{At}y_{0} = Ae^{At}y_{0} = Ay. label{4.14}]

Par conséquent, le problème général de résolution (4.10) est équivalent au calcul de (e^{At}), et nous tournons maintenant notre attention vers cette tâche.

Tout d'abord, supposons que (A) est une matrice diagonale, disons

[A = egin{pmatrix} {lambda_{1}}&{0}&{cdots}&{0} {0}&{lambda_{2}}&{cdots}&{0 } {0}&{0}&{cdots}&{0} {0}&{0}&{cdots}&{lambda_{n}} end{pmatrix} label{4.15 }]

Ensuite, il est facile de voir en remplaçant UNE dans la série exponentielle (4.12) qui :

[e^{At} = egin{pmatrix} {e^{lambda_{1}t}}&{0}&{cdots}&{0} {0}&{e^{lambda_ {2}t}}&{cdots}&{0} {0}&{0}&{cdots}&{0} {0}&{0}&{cdots}&{e ^{lambda}t} end{pmatrix} label{4.16}]

Par conséquent, notre stratégie sera de transformer les coordonnées de sorte que dans les nouvelles coordonnées A devienne diagonale (ou le plus "proche que possible" de la diagonale, ce que nous expliquerons bientôt). Alors (e^{At}) sera facilement calculable dans ces coordonnées. Une fois cela accompli, nous utilisons l'inverse de la transformation pour retransformer la solution dans le système de coordonnées d'origine.

Maintenant, précisons ces idées. Nous laissons

[y = Tu, u in mathbb{R}^n, y in mathbb{R}^n, label{4.17}]

où T est une matrice (n imes n) dont les propriétés précises seront développées dans la suite.

C'est une approche typique dans les EDO. Nous proposons une transformation générale des coordonnées de l'ODE, puis nous la construisons d'une manière qui donne les propriétés de l'ODE que nous désirons. La substitution de (4.17) en (4.10) donne :

[dot{y} = Tdot{u} = Ay = ATu, label{4.18}]

T sera construit de manière à le rendre inversible, de sorte que nous avons :

[dot{u} = T^{-1}ATu, u(0) = T^{-1}y(0). label{4.19}]

Pour simplifier la notation on laisse :

[T = T^{-1}AT, label{4.20}]

ou alors

[A = T^{-1}Lambda T. label{4.21}]

En substituant (4.21) dans la série l'exponentielle matricielle (4.12) donne :

(e^{At} = e^{TLambda T^{-1}t}),

[= mathbb{1}+ TLambda T^{-1}t+frac{1}{2!}(TLambda T^{-1})^{2}t^2+cdots+ frac{1}{n!}(TLambda T^{-1})^{n}t^n+cdots label{4.22}]

Notez maintenant que pour tout entier positif n nous avons :

((TLambda T^{-1})^{n} = underbrace{(TLambda T^{-1})(TLambda T^{-1})cdots (TLambda T ^{-1})(TLambda T^{-1})}_{n facteurs})

[= TLambda^{n}T^{-1} label{4.23}]

Substituer ceci dans (4.22) donne :

(e^{At} = sum_{n=0}^{infty} frac{1}{n!}(TLambda T^{-1})^{n}t^n),

(= T(sum_{n=0}^{infty} frac{1}{n!}Lambda^{n}t^n)T^{-1}),

[= Te^{Lambda t}T^{-1} label{4.24}]

ou alors

[e^{At} = Te^{Lambda t}T^{-1} label{4.25}]

Nous arrivons maintenant à notre résultat principal. Si T est construit de telle sorte que

[Lambda = T^{-1}AT label{4.26}]

est diagonale, alors il résulte de (4.16) et (4.25) que (e^{At}) peut toujours être calculé. Ainsi, le problème ODE de résolution (4.10) devient un problème d'algèbre linéaire. Mais une matrice (n imes n) générale A peut-elle toujours être diagonalisée ? Si vous avez suivi un cours d'algèbre linéaire, vous savez que la réponse à cette question est « non ». Il y a une théorie du (réel) qui s'appliquera ici. Cependant, cela nous entraînerait dans une trop grande diversion pour ce cours. Au lieu de cela, nous considérerons les trois cas standard pour les matrices (2 imes 2). Cela suffira pour introduire les idées principales sans s'enliser dans l'algèbre linéaire. Néanmoins, il ne peut pas être totalement évité. Vous devrez être capable de calculer les valeurs propres et les vecteurs propres des matrices (2 imes 2), et comprendre leur signification.

Les trois cas de matrices (2 imes 2) que nous allons considérer sont caractérisés par leurs valeurs propres :

  • deux valeurs propres réelles, diagonalisables A,
  • deux valeurs propres identiques, A non diagonalisable,
  • une paire conjuguée complexe de valeurs propres.

Dans le tableau ci-dessous, nous résumons la forme sous laquelle ces matrices peuvent être transformées (appelée la de A) et l'exponentielle résultante de cette forme canonique.

Une fois la transformation en (Lambda) effectuée, nous utiliserons ces résultats pour en déduire (e^{Lambda}).


4 : Comportement à proximité de trajectoires - Linéarisation - Mathématiques

Math 321 : Équations différentielles Printemps 2005

Dr James R. Hughes
Département des sciences mathématiques du Collège Elizabethtown

6.2 Systèmes linéaires et presque linéaires

  • Examen de la classification de type et de stabilité pour les systèmes linéaires à coefficient constant 2 par 2
  • Changement de coordonnées pour déplacer un point critique vers l'origine
  • Systèmes presque linéaires
    • Définitions : point critique isolé, système linéarisé, système presque linéaire
    • Théorème 2 : Stabilité des systèmes presque linéaires
    • Linéarisation près d'un point critique isolé (exemple)

    Examen de la classification de type et de stabilité pour les systèmes linéaires à coefficient constant 2 par 2

    x'(t) = a x + b y
    y'(t) = c x + d y

    (où a, b, c et d sont des constantes) a un seul point critique à (0,0) et le type et la stabilité de (0,0) peuvent être complètement déterminés à partir des valeurs propres de la matrice de coefficients. Spécifiquement:

    S'il existe deux valeurs propres réelles de signe opposé, l'origine est un point selle (et donc instable).

    S'il y a deux valeurs propres réelles positives, l'origine est un nœud instable.

    S'il y a deux valeurs propres réelles négatives, l'origine est un nœud asymptotiquement stable.

    Si les valeurs propres sont complexes de partie réelle négative, l'origine est un point spiral asymptotiquement stable.

    Si les valeurs propres sont complexes à partie réelle positive, l'origine est un point spiral instable.

    Et enfin, si les valeurs propres sont complexes à partie réelle nulle, l'origine est un centre (et donc stable).

    Le texte fait une distinction entre correct et non conforme nœuds. En bref, un nœud est correct si toutes les trajectoires s'approchent ou s'éloignent de l'origine le long de lignes tangentes différentes, et non conforme si toutes les trajectoires s'approchent ou s'éloignent le long de la même tangente. La plupart des nœuds que nous verrons sont impropres. Un nœud propre n'apparaît que dans le cas d'une valeur propre complète de multiplicité 2 (voir la figure 6.2.3 dans votre texte).

    Changement de coordonnées pour déplacer un point critique vers l'origine

    Si un système autonome a un point critique en (x0, y0), alors le changement de variables u = x - x0, v = y - y0 aboutit à un système équivalent à l'original, mais traduit de sorte que le point critique à (x0, y0) est déplacé vers (0, 0).

    x'(t) = 3x - 4y - 2
    y'(t) = 2x + y - 5

    a un seul point critique en (2, 1). En fixant u = x - 2 et v = y - 1, on a
    u' = x' , v' = y', x = u + 2, et y = v + 1, et le système devient

    u'(t) = 3(u + 2) - 4(v + 1) - 2
    v'(t) = 2(u + 2) + (v + 1) - 5

    qui a un seul point critique à (0, 0). Le type et la stabilité du point critique d'origine peuvent maintenant être déterminés à partir de ceux de (0,0) dans le nouveau système. Le polynôme caractéristique du nouveau système est
    (3 - r)(1 - r) + 8 = r 2 - 4r + 11, qui a des racines complexes avec une partie réelle positive, donc le point critique est un point spiral instable. Cela peut être confirmé en examinant un portrait de phase généré par Maple pour le système d'origine :


    Nous allons maintenant considérer le système autonome plus général (*) ci-dessous.

    Un point critique (x0, y0) de (*) est appelé isolé si un rectangle ouvert le contenant ne contient aucun autre point critique.

    Si (x0, y0) est un point critique isolé, alors, comme ci-dessus, le changement de variables u = x - x0, v = y - y0 convertit (*) en un système dans lequel le point critique en (x0, y0) devient un point critique à (0, 0). De plus, si le système transformé peut s'écrire sous la forme

    u'(t) = a u + b v + f(u, v)
    v'(t) = c u + d v + g(u, v)

    où f(u, v) et g(u, v) ont la propriété que les limites lorsque (u, v) s'approchent de (0, 0) des deux

    f(u, v)/sqrt(u 2 + v 2 ) et g(u, v)/sqrt(u 2 + v 2 )

    sont 0 (comme cela se produirait si f et g n'étaient constitués que de puissances de u et v supérieures à 1),
    alors le système est appelé presque linéaire et le système linéaire

    u'(t) = a u + b v
    v'(t) = c u + d v

    est appelé le linéarisation de (*). Dans la plupart des conditions, le comportement qualitatif des trajectoires proches de (0,0) dans le système linéarisé ressemble à celui des trajectoires proches de (x0, y0) dans le système d'origine. Concrètement, nous avons

    Théorème 2 (Stabilité des systèmes presque linéaires)

    (a) Si la matrice de coefficients de la linéarisation de (*) a une seule valeur propre réelle r de multiplicité 2, alors le point critique (x0, y0) est soit un nœud, soit un point en spirale, et est asymptotiquement stable si r négatif, et instable si r est positif.

    (b) Si la matrice de coefficients de la linéarisation de (*) a des valeurs propres imaginaires pures, alors le point critique (x0, y0) est soit un centre, soit un point en spirale, et peut être asymptotiquement stable, stable ou instable.

    (c) Dans tous les autres cas, le type et la stabilité du point critique (x0, y0) de (*) sont les mêmes que ceux de (0,0) dans la linéarisation de (*).

    x'(t) = x 2 + 3xy + 2y 2
    y'(t) = 4 - x 2 .

    Nous trouvons d'abord les points critiques. Le réglage y' = 0 nous indique que x doit être 2 ou -2. La substitution de x = 2 dans la première équation, et la définition de x' = 0, nous indiquent que y doit être -1 ou -2, donc deux points critiques sont (2, -1) et (2, -2). De même, la substitution de x = -2 dans la première équation nous permet d'avoir deux autres points critiques, (-2, 1) et (-2, 2).

    Nous linéarisons d'abord en (2, -1). On pose u = x - 2 et v = y + 1, donc u' = x', v' = y', x = u + 2, et y = v - 1, et le système est transformé en

    u'(t) = (u + 2) 2 + 3(u + 2)(v - 1) + 2(v - 1) 2
    v'(t) = 4 - (u + 2) 2 .

    En multipliant les côtés droits, on obtient

    u'(t) = u + 2v + u 2 + 3uv + 2v 2
    v'(t) = -4u - u 2

    La matrice de coefficients du système linéarisé a des valeurs propres complexes avec une partie réelle positive, donc (2, -1) est un point de spirale instable.

    Ensuite, nous linéarisons en (2, -2). En fixant u = x - 2 et v = y + 2, le système est transformé en

    u'(t) = (u + 2) 2 + 3(u + 2)(v - 2) + 2(v - 2) 2
    v'(t) = 4 - (u + 2) 2

    u'(t) = -2u - 2v + u 2 + 3uv + 2v 2
    v'(t) = -4u - u 2 ,

    Les valeurs propres de la matrice de coefficients sont 2 et -4 dans ce cas, donc le point critique (2, -2) du système d'origine est un point de selle instable.

    La linéarisation à (-2, 1) donne

    u'(t) = (u - 2) 2 + 3(u - 2)(v + 1) + 2(v + 1) 2
    v'(t) = 4 - (u - 2) 2

    u'(t) = -u - 2v + u 2 + 3uv + 2v 2
    v'(t) = 4u - u 2 ,

    Les valeurs propres de la matrice de coefficients sont complexes avec une partie réelle négative dans ce cas, donc le point critique (-2, 1) du système d'origine est un point spiral asymptotiquement stable.

    Enfin, nous linéarisons en (-2, 2). Réglage u = x + 2 et v = y - 2

    u'(t) = (u - 2) 2 + 3(u - 2)(v + 2) + 2(v + 2) 2
    v'(t) = 4 - (u - 2) 2

    u'(t) = 2u + 2v + u 2 + 3uv + 2v 2
    v'(t) = 4u - u 2 ,

    La matrice de coefficients a des valeurs propres 4 et -2, donc encore une fois le point critique (-2,2) est un point selle instable.


    Notes sur Diffy Qs : équations différentielles pour les ingénieurs

    Remarque : 1 cours magistral, §6.1–§6.2 dans [EP], §9.2–§9.3 dans [BD]

    À l'exception de quelques brefs détours au chapitre 1, nous avons considéré principalement des équations linéaires. Les équations linéaires suffisent dans de nombreuses applications, mais en réalité la plupart des phénomènes nécessitent des équations non linéaires. Les équations non linéaires, cependant, sont notoirement plus difficiles à comprendre que les équations linéaires, et de nombreux nouveaux phénomènes étranges apparaissent lorsque nous permettons à nos équations d'être non linéaires.

    Ne vous inquiétez pas, nous n'avons pas perdu tout ce temps à étudier des équations linéaires. Les équations non linéaires peuvent souvent être approchées par des équations linéaires si nous n'avons besoin d'une solution « localement », par exemple, que pour une courte période de temps, ou uniquement pour certains paramètres. Comprendre les équations linéaires peut également nous donner une compréhension qualitative d'un problème non linéaire plus général. L'idée est similaire à ce que vous avez fait en calcul en essayant d'approcher une fonction par une droite avec la bonne pente.

    Dans la section 2.4, nous avons examiné le pendule de longueur (L ext<.>) Le but était de résoudre pour l'angle ( heta(t)) en fonction du temps (t ext<. >) L'équation pour la configuration est l'équation non linéaire

    Au lieu de résoudre cette équation, nous avons résolu l'équation linéaire plutôt simple

    Bien que la solution de l'équation linéaire ne soit pas exactement ce que nous recherchions, elle est plutôt proche de l'originale, tant que l'angle ( heta) est petit et que la période de temps impliquée est courte.

    Vous pourriez demander : pourquoi ne résolvons-nous pas simplement le problème non linéaire ? Eh bien, cela peut être très difficile, peu pratique ou impossible à résoudre analytiquement, selon l'équation en question. Nous pouvons même ne pas être intéressés par la solution réelle, nous pourrions seulement être intéressés par une idée qualitative de ce que fait la solution. Par exemple, que se passe-t-il lorsque le temps tend vers l'infini ?

    Sous-section 8.1.1 Systèmes autonomes et analyse du plan de phase

    Nous limitons notre attention à un système autonome à deux dimensions

    où (f(x,y)) et (g(x,y)) sont des fonctions de deux variables, et les dérivées sont prises par rapport au temps (t ext<.>) Les solutions sont des fonctions (x(t)) et (y(t)) tels que

    La façon dont nous analyserons le système est très similaire à la section 1.6, où nous avons étudié une seule équation autonome. Les idées en deux dimensions sont les mêmes, mais le comportement peut être beaucoup plus compliqué.

    Il peut être préférable de considérer le système d'équations comme l'équation vectorielle unique

    Comme dans la section 3.1, nous dessinons le portrait de phase (ou alors diagramme de phase), où chaque point ((x,y)) correspond à un état spécifique du système. Nous dessinons le champ vectoriel donné en chaque point ((x,y)) par le vecteur (left[ egin f(x,y) g(x,y) end ight] ext<.>) Et comme avant si on trouve des solutions, on trace les trajectoires en traçant tous les points (igl(x(t),y(t)igr)) pour une certaine plage de (t exte<.>)

    Exemple 8.1.1.

    Considérons l'équation du second ordre (x''=-x+x^2 ext<.>) Écrivez cette équation comme un système non linéaire du premier ordre

    Le portrait de phase avec quelques trajectoires est tracé à la figure 8.1.

    Graphique 8.1. Portrait de phase avec quelques trajectoires de (x' = y ext<,>) (y' = -x+x^2 ext<.>)

    D'après le portrait de phase, il devrait être clair que même ce système simple a un comportement assez compliqué. Certaines trajectoires continuent d'osciller autour de l'origine, d'autres partent vers l'infini. Nous reviendrons souvent sur cet exemple et l'analyserons complètement dans cette (et la suivante) section.

    Si nous zoomons sur le diagramme près d'un point où (left[ egin f(x,y) g(x,y) end ight]) n'est pas zéro, alors à proximité, les flèches pointent généralement dans essentiellement la même direction et ont essentiellement la même amplitude. En d'autres termes, le comportement n'est pas si intéressant près d'un tel point. Nous supposons bien sûr que (f(x,y)) et (g(x,y)) sont continus.

    Concentrons-nous sur les points du diagramme de phases ci-dessus où les trajectoires semblent commencer, se terminer ou faire le tour. Nous voyons deux de ces points : ((0,0)) et ((1,0) ext<.>) Les trajectoires semblent faire le tour du point ((0,0) ext<,> ) et ils semblent entrer ou sortir du point ((1,0) ext<.>) Ces points sont précisément les points où les dérivées de (x) et (y) sont nuls. Définissons le points critiques comme les points ((x,y)) tels que

    En d'autres termes, ce sont les points où (f(x,y)=0) et (g(x,y)=0 ext<.>)

    Les points critiques sont là où le comportement du système est en quelque sorte le plus compliqué. Si (gauche[ egin f(x,y) g(x,y) end ight]) vaut zéro, alors à proximité, le vecteur peut pointer dans n'importe quelle direction. De plus, les trajectoires vont soit vers, loin ou autour de ces points, donc si nous recherchons un comportement qualitatif à long terme du système, nous devrions regarder ce qui se passe près des points critiques.

    Les points critiques sont aussi parfois appelés équilibres, puisque nous avons ce qu'on appelle solutions d'équilibre aux points critiques. Si ((x_0,y_0)) est un point critique, alors nous avons les solutions

    Dans l'exemple 8.1.1, il existe deux solutions d'équilibre :

    Comparez cette discussion sur les équilibres à celle de la section 1.6. Le concept sous-jacent est exactement le même.

    Sous-section 8.1.2 Linéarisation

    Dans la section 3.5, nous avons étudié le comportement d'un système linéaire homogène de deux équations près d'un point critique. Pour un système linéaire de deux variables donné par une matrice inversible, le seul point critique est l'origine ((0,0) ext<.>) Mettons à profit la compréhension que nous avons acquise dans cette section pour comprendre ce qui se passe près des points critiques des systèmes non linéaires.

    En calcul, nous avons appris à estimer une fonction en prenant sa dérivée et en la linéarisant. Nous travaillons de la même manière avec des systèmes non linéaires d'EDO. Supposons que ((x_0,y_0)) soit un point critique. Remplacez d'abord les variables par ((u,v) ext<,>) de sorte que ((u,v)=(0,0)) corresponde à ((x_0,y_0) ext<.> ) C'est-à-dire,

    Ensuite, nous devons trouver la dérivée. Dans le calcul multivariable, vous avez peut-être vu que la version à plusieurs variables de la dérivée est la matrice jacobienne 1 . La matrice Jacobienne de la fonction à valeur vectorielle (left[ egin f(x,y) g(x,y) end ight]) à ((x_0,y_0)) est

    Cette matrice donne la meilleure approximation linéaire car (u) et (v) (et donc (x) et (y)) varient. Nous définissons le linéarisation de l'équation (8.1) comme système linéaire

    Exemple 8.1.2.

    Restons avec les mêmes équations que l'exemple 8.1.1 : (x' = y ext<,>) (y' = -x+x^2 ext<.>) Il y a deux points critiques, ((0,0)) et ((1,0) ext<.>) La matrice Jacobienne en tout point est

    Donc en ((0,0) ext<,>) on a (u=x) et (v=y ext<,>) et la linéarisation est

    Au point ((1,0) ext<,>) nous avons (u=x-1) et (v=y ext<,>) et la linéarisation est

    Les diagrammes de phase des deux linéarisations au point ((0,0)) et ((1,0)) sont donnés sur la figure 8.2. Notez que les variables sont maintenant (u) et (v ext<.>) Comparez la figure 8.2 avec la figure 8.1, et regardez surtout le comportement près des points critiques.

    Graphique 8.2. Diagramme de phases avec quelques trajectoires de linéarisations aux points critiques ((0,0)) (gauche) et ((1,0)) (droite) de (x' = y ext<,>) (y' = -x+x^2 exte<.>)


    Portraits de phase des systèmes non linéaires

    Considérez un , éventuellement non linéaire, système autonome , (autonome signifie que la variable indépendante , considéré comme représentant le temps, ne se produit pas du côté droit des équations). Tout comme nous l'avons fait pour les systèmes linéaires, nous voulons regarder les trajectoires du système. Comme dans la plupart des cas il est impossible de résoudre exactement ces systèmes, nous nous concentrerons sur les aspects qualitatifs, et en particulier sur la manière de tracer le portrait de phase à la main.
    • Les trajectoires suivent le champ de direction. Le vecteur vitesse pour une solution en un point dans l'avion est . La direction de la trajectoire est la direction de ce vecteur.
    • Les courbes et sont les isoclines sur lesquelles la direction d'une trajectoire est respectivement verticale et horizontale. Ces isoclines divisent le plan en régions. Dans chaque région, les signes de et ne changez pas, donc par ex. si les deux sont positifs, la direction est toujours vers le haut et vers la droite. Généralement, lorsque vous traversez une isocline, une composante de la vitesse change de signe. Ainsi, en connaissant la direction en un point, vous pouvez déterminer les directions dans toutes les régions.
    • Une isocline (ou une partie de celle-ci) qui est une ligne verticale est une trajectoire (ou peut-être plusieurs d'entre elles). De même, une isocline (ou une partie de celle-ci) qui est une ligne horizontale est une trajectoire (ou plusieurs).
    • Les trajectoires ne se rencontrent ni ne s'arrêtent, sauf qu'à la limite comme ou alors ils peuvent s'approcher d'un point d'équilibre.
    • Les points d'équilibre (également appelés points critiques ou stationnaires) sont les points où les deux et . Ils sont donc aux intersections de ces isoclines.
    • Le comportement à proximité des points d'équilibre est important. Les trajectoires proches d'un point d'équilibre sont très bien approchées par celles de la linéarisation du système en ce point, et le point critique peut être classé grâce à cette linéarisation (à deux exceptions près que nous verrons).
    • Une séparatrice (pluriel séparatrices ) est une trajectoire séparant deux régions dans lesquelles le comportement des solutions comme ou alors est différent. Les trajectoires entrant et sortant d'un point selle sont très souvent séparatrices. Ceux-ci devraient être parmi les premières trajectoires que vous esquissez.
    • Si le système a une certaine symétrie, cela peut aider. Des exemples de symétrie incluent l'échange et , ou changer le signe de et/ou .


    L'approximation linéaire d'une fonction de deux variables près d'un point est

    Nous voulons le faire en particulier à un point d'équilibre, où . Il est utile de faire un changement de variables , , alors , correspond au point d'équilibre. Alors la linéarisation de notre système à ce point d'équilibre est la système homogène linéaire à coefficient constant

    est appelé le Jacobien du champ de vecteurs .

    Exemple : Considérez le système

    Pour un point d'équilibre, nous avons besoin (c'est à dire. ou alors ) et (c'est à dire. ou alors ). Il y a donc deux points d'équilibre : et .

    Avant de classer les points d'équilibre, c'est une bonne idée de tracer les isoclines pour et . Je tracerai le premier en bleu et le second en vert. J'indique avec des flèches le champ de direction dans chaque région et sur les isoclines. Les intersections des isoclines bleues et vertes sont les points d'équilibre.

    La matrice Jacobienne est . Au point d'équilibre c'est . Qui a une double valeur propre , et deux vecteurs propres linéairement indépendants. Donc le point d'équilibre est un nœud singulier, et est un attracteur.

    Au deuxième point d'équilibre la matrice Jacobienne est . Celui-ci a les valeurs propres 1 et . Donc le point d'équilibre est une selle. Les vecteurs propres sont pour 1 et pour . L'image ci-dessous montre le plan de phase avec certaines parties des trajectoires à proximité des deux points d'équilibre. Notez que les directions de ces trajectoires concordent avec les flèches du champ de direction de l'image précédente.

    Esquissons maintenant quelques trajectoires. Il y a des trajectoires sur le et axes (vers l'intérieur jusqu'au point d'équilibre à l'origine), car lorsque et lorsque . Ensuite, c'est une bonne idée d'esquisser les trajectoires sortant et entrant au point de selle. Par exemple, on arrive au point de selle par le bas et vers la droite. On remonte le temps en suivant les flèches en arrière. Ces flèches pointent vers le haut et vers la gauche dans la région , . La trajectoire doit courber pour éviter la trajectoire sur le positif axe. Il est vraisemblablement asymptotique à cet axe. La trajectoire sortant du point de selle vers le bas et vers la gauche doit continuer vers le bas et vers la gauche jusqu'à se terminer au point d'équilibre .

    Enfin, nous dessinons d'autres trajectoires, dont au moins une dans chaque région. Notez que les trajectoires qui entrent dans le point d'équilibre à peut le faire à n'importe quel angle.

    Notre image est symétrique par rapport à la ligne , en raison du fait que le système d'équations reste le même si vous intervertissez et .

    Les trajectoires entrant et sortant du point selle sont séparatrices. Toutes les trajectoires en dessous et à gauche des deux trajectoires entrant dans la selle vont au point d'équilibre comme , tandis que ceux au-dessus et à droite s'éloignent à l'infini asymptotique à la ligne . En bas et à droite des trajectoires de sortie de selle, tout rentre de l'infini asymptotique au positif axe (comme ), tandis qu'au-dessus et à gauche de ceux-ci tout vient de l'infini asymptotique au positif axe.


    Comme je l'ai mentionné, il existe deux exceptions à la règle selon laquelle le portrait de phase près d'un point d'équilibre peut être classé par la linéarisation à ce point d'équilibre. La première est où 0 est une valeur propre de la linéarisation (nous n'avons même pas regardé le système linéaire dans ce cas !). La deuxième exception est lorsque la linéarisation est un centre. Le système linéaire a des solutions périodiques, correspondant à des trajectoires qui sont des courbes fermées (ellipses). Imaginez commencer à un moment donné et suivre le champ de direction. Après avoir fait tout le tour du point d'équilibre, dans le système linéaire, vous revenez exactement au point de départ. C'est une question très délicate, et tout effet non linéaire, même minime, pourrait la gâcher. Si dans le système non linéaire vous revenez un peu plus loin du point d'équilibre que là où vous avez commencé, alors votre trajectoire ne peut pas être une courbe fermée. La prochaine fois, vous serez encore plus loin. La trajectoire s'éloignera du point d'équilibre en spirale. Si toutes les trajectoires proches du point d'équilibre sont comme ceci, le point d'équilibre est une spirale instable. En revanche, si après un tour du point d'équilibre vous êtes un peu plus près que votre point de départ, votre trajectoire s'incurve vers l'intérieur. Si toutes les trajectoires proches du point d'équilibre sont comme ceci, le point d'équilibre est une spirale stable (et un attracteur). Voici des photos de ces deux possibilités. La troisième possibilité, bien sûr, est que vous reveniez exactement au point où vous avez commencé, et c'est vraiment un centre.

    Une façon de montrer qu'un centre du système linéarisé est toujours un centre dans le système non linéaire est de trouver une équation pour les trajectoires. S'il existe une telle équation (implicite) est une fonction lisse, non constante dans aucune région, et une constante arbitraire, alors les trajectoires, qui sont des courbes de niveau de cette fonction, ne peuvent pas être des spirales mais peuvent être des courbes fermées. Cela se produit à la fois dans les systèmes prédateur-proie et pendule.


    3 réponses 3

    En général, vous pouvez linéariser autour de n'importe quel solution connue. L'idée est qu'une fois qu'une solution $ heta_0(t)$ est connue, les solutions voisines $ heta$ suivent approximativement une équation linéaire : à savoir, écrire $ heta_0(t) = heta_0(t) + h(t)$ on obtient $ I heta''+Mglsin heta approx (I heta_0''+Mglsin heta_0) + Ih''+(Mglcos heta_0 )h $ ce qui conduit à une équation linéaire approximative $ Ih''+(Mglcos heta_0 )h =0$ car $(I heta_0''+Mglsin heta_0)=u$.

    Le hic, c'est : connaissez-vous $ heta_0$ pour commencer ? Une solution d'équilibre est facile à trouver. Trouver une solution générique. eh bien, c'est juste le problème d'origine.

    But you will occasionally see linearization along a non-constant periodic orbit called a limit cycle or even an arbitrary trajectory. This is generally referred to as tracking-control.


    Glossaire

    • g(X) is an analytic function at the origin (or, more precise, admits the second order Taylor's approximation)
    • as X &rarr 0,

    Example 1: Not almost linear system

    To analyze the trajectories of Eq.eqref, it is convenient to use polar coordinates

    If there exists a &delta: 0 < &delta &le r₀, such that, for any solution path &phi = &langle &phi₁, &phi₂ &rangle of nonlinear system eqref that has at least one point inside the circle 0 < r < &delta, the solution exists over a t half line, and if


    1 Answer 1

    OK, here's a quick rundown on how I would do this hopefully our OP Mirjan Pecenko can check his own work against what I do here.

    $dot x = xy + 1 = dot y = xy + x = 0 ag 3$

    $xy = -x Longrightarrow -x + 1 = 0 Longrightarrow x = 1, ag 4$

    $xy + 1 = 0 Longrightarrow xy = -1 Longrightarrow y = (1)y = -1 ag 4$

    so it appears the only critical point is at $(1, -1)$. The Jacobean of the vector field

    $vec V(x, y) = egin xy + 1 x + xy end ag 5$

    $J_V(1, -1) = egin (partial(xy + 1)/partial x & (partial(xy + 1)/partial y (partial(xy + x)/partial x & (partial(xy + x)/partial y end_<(1, -1)>$ $= egin y & x y + 1 & x end_ <(1, -1)>= egin -1 & 1 0 & 1 end ag 6$

    it is now obvious that the eigenvalues of $J_V(1, -1)$ are $pm 1$ therefore this point is a saddle the eigenvectors at $(1, -1)$ are $(1, 0)^T$ for $-1$ and $(1/2, 1)^T$ for $1$ it is now easy to sketch a phase portrait for this system, a task I leave to my readers.

    It should be noted that when sketching a phase portrait, it is often helpful to find those curves in $Bbb R^2$ where $dot x = 0$ and/or $dot y = 0$. These curves are useful guides to finding the geometry of the solutions, since they show us where the tangent lines to the integral curves or vertical or horizontal, respectively. When combined with the shapes of the solutions near the critical point provided by the above analysis, we can get a pretty good idea of how the flow will appear. As with any hand-done graphical analysis, care must be taken to ensure that we draw accurately enough to capture only correct features of the trajectories.

    Note Added in Edit, Thursday 14 June 2018 12:35 PM PST: It appears that the notion of isoclines, which proves to be most convenient as a guide to sketching phase portraits, may be generalized in a way which allows the extraction of more information about the integral curves and/or overall shape of the flow of a given vector field. In this problem, isoclines are only exploited in a rather peripheral way, since they are merely mentioned as a sort of after-thought in the comments. Nevertheless, they may be used much more extensively. Indeed, rather than restricting the use of the isocline method to determining the curves on which $dot x = 0$ and/or $dot y = 0$, we may if we so choose deploy it in an attempt to find just where $dot x, dot y$ take on any of their possible values. One technique which can help effect this aim is to use the gradient of the component functions, in this case $xy + 1$ and $xy + x$, to guide us towards regions of greater or lesser component magnitude for example, since

    which indicates the direction in which $dot x$ increases, so that we may, for example, find the directions in which points on the $dot x = 0$ isocline must be moved to make $dot x$ larger. By systematic application of such methodology, quite detailed phase portraits may be had. Unfortunately, at present I lack both the graphics tools and the time to illustrate what I am saying passant par pictorial means. End of Note.


    Using basic algebra, we find that the critical points of this system are (0,0), (1,0.25), (1.125,0). These points are marked with red stars on the direction field in part (a).

    Critical Point (0,0)

    Using formula (13) from Section 9.3, we find that the linear system corresponding to the almost linear system at (0,0) is given by the matrix equation

    The eigenvalues and eigenvectors for this linear system are:

    The eigenvalues are real and have opposite signs, so the origin is a saddle point and is thus unstable.

    Critical Point (1.125, 0)

    The almost linear system at (1.125, 0) is given by the matrix equation

    The corresponding eigenvalues for this linear system are:

    The eigenvalues are real and have opposite signs, so this point is also a saddle point and is unstable.

    Critical Point (1, 0.25)

    Again using formula (13) from Section 9.3, we see that linear system that approximates the almost linear system at the critical point (1, 0.25) is:

    This linear system has eigenvalues and eigenvectors:

    The eigenvalues for this linear system are both negative and are unequal, so the point (1, 0.25) is an asymptotically stable node.


    Glossaire

    Some situations require more than one differential equation to model a particular situation. We might use a system of differential equations to model two interacting species, say where one species preys on the other. For example, we can model how the population of Canadian lynx (lynx Canadians) interacts with a the population of snowshoe hare (lepus americanis) (see https://www.youtube.com/watch?v=ZWucOrSOdCs).

    A good historical data are known for the populations of the lynx and snowshoe hare from the Hudson Bay Company. This large Canadian retail company, which owns and operates a large number of retail stores in North America and Europe, including Galeria Kaufhof, Hudson's Bay, Home Outfitters, Lord & Taylor, and Saks Fifth Avenue, was originally founded in 1670 as a fur trading company. The Hudson Bay Company kept accurate records on the number of lynx pelts that were bought from trappers from 1821 to 1940. The company noticed that the number of pelts varied from year to year and that the number of lynx pelts reached a peak about every ten years. The ten year cycle for lynx can be best understood using a system of differential equations.

    The primary prey for the Canadian lynx is the snowshoe hare. We will denote the population of hares by H(t) and the population of lynx by L(t), où t is the time measured in years. We will make the following assumptions for our predator-prey model.

      If no lynx are present, we will assume that the hares reproduce at a rate proportional to their population and are not affected by overcrowding. That is, the hare population will grow exponentially,

    The Lotka--Volterra system of equations is an example of a Kolmogorov model, which is a more general framework that can model the dynamics of ecological systems with predator-prey interactions, competition, disease, and mutualism. The equations which model the struggle for existence of two species (prey and predators) bear the name of two scientists: Lotka (1880--1949) and Volterra (1860--1940). They lived in different countries, had distinct professional and life trajectories, but they are linked together by their interest and results in mathematical modeling.

    The predator–prey model was initially proposed by Alfred J. Lotka in the theory of autocatalytic chemical reactions in 1910. Lotka was born in Lemberg, Austria-Hungary, but his parents immigrated to the US. In 1925, he utilized the equations to analyze predator-prey interactions. Lotka published almost a hundred articles on various themes in chemistry, physics, epidemiology or biology, about half of them being devoted to population issues. He also wrote six books.

    The same set of equations was published in 1926 by Vito Volterra, a mathematician and physicist, who had become interested in mathematical biology because of the impact by the marine biologist Umberto D'Ancona (1896--1964). Vito Volterra was born in Ancona, then part of the Papal States, into a very poor Jewish family. He attended the University of Pisa, where he became professor of rational mechanics in 1883. His most famous work was done on integral equations. In 1892, he became professor of mechanics at the University of Turin and then, in 1900, professor of mathematical physics at the University of Rome La Sapienza. His daughter Luisa married Umberto D’Ancona.

    The predator-prey system of equations was later extended by many researchers, including C. S. Holling, Arditi--Ginzburg model, Rosenzweig-McArthur model, and some others.

    The critical points of the Lotka--Volterra system of equations are the solutions of the algebraic equations

    We may try to find the general solution of the Lotka--Volterra system of equations. From both equations, we get

    Notice that the predator population, L, begins to grow and reaches a peak after the prey population, H reaches its peak. As the prey population declines, the predator population also declines. Once the predator population is smaller, the prey population has a chance to recover that the cycle begins again.

    we can place two plots sude by side:

    1. ( displaystyle 0 < a_1 delta_1 < K left( m_1 - delta_1 ight) ) and
    2. ( displaystyle a_2 delta_2 > K left( m_2 - delta_2 ight) ) or ( b_2 le 1 . )

    Mathematical analysis of the Beddington--DeAngelis system shows that there exist two equilibria (0,0) and ( left( frac , frac ight) = (4,1) , ) being globally stable in the interior of the first quadrant. The eigenvalues of the Jacobian matrix at the origin are ( lambda_1 =1 quadmboxquad lambda_2 =5 , ) and the eigenvalues of the Jacobian at the point (4,1) are given by ( - frac<1> <12>pm <f j>, frac> <12>. )

    1. Barabas, Gyorgy, D'Andrea, Rafael, and Stump, Simon Maccracken, Chesson's coexistence theory, Ecological Monographs, 2018, https://doi.org/10.1002/ecm.1302
    2. Batiha, K., Numerical Solutions of the Multispecies Predator-Prey Model by Variational Iteration Method, Journal of Computer Science, 2007, Vol. 3 (7): 523-527, 2007
    3. Bayliss, A., Nepomnyashchy, A.A., Volpert, V.A., Mathematical modeling of cyclic population dynamics, Physica D: Nonlinear Phenomena, 2019, Volume 394, doi: 10.1016/j.physd.2019.01.010
    4. Dellal, M., Lakrib, M., Sari, T., The operating diagram of a model of two competitors in a chemostat with an external inhibitor, Mathematical Biosciences, 302, No 8, 2018, 27--45.
    5. Dimitrov, D.T. and Kojouharov, H.V., Complete mathematical analysis of predator–prey models with linear prey growth and Beddington–DeAngelis functional response, Applied Mathematics and Computation, 162, (2), 523--538, 2005.
    6. Goel, N.S., Maitra, S.C., and Montroll, E., On the Volterra and other nonlinear models of interacting populations, Reviews of Modern Physics, 1971, Vol. 43, pp. 231-- https://doi.org/10.1103/RevModPhys.43.231
    7. Hsu, S.B., Hubbell, S.P., and Paul Waltman, Competing predators, SIAM Journal on Applied Mathematics, 35, No 4, 1978, 617--625.
    8. Hsu, S.B., Hubbell, S.P., and Paul Waltman, A mathematical theory for single-nutrient competition in continuous cultures of micro-organisms, SIAM Journal on Applied Mathematics, Vol 32, No 2, 1977, 366--383.
    9. Hsu, S.B., Hubbell, S.P., and Paul Waltman, A contribution to the theory of competing predators, Ecological Monographs, Vol 48, No 3, 1978, 337--349.
    10. May, R.M., Leonard, W.J., Nonlinear Aspects of Competition Between Three Species, SIAM Journal on Applied Mathematics, Vol. 29, Issue 2, pp. https://doi.org/10.1137/0129022
    11. Molla, H., Rahman, M.S., Sarwardi, S., Dynamics of a Predator–Prey Model with Holling Type II Functional Response Incorporating a Prey Refuge Depending on Both the Species, International Journal of Nonlinear Sciences and Numerical Simulation, 2018, Vol. 20, No. 1, https://doi.org/10.1515/ijnsns-2017-0224
    12. Olek. S. An accurate solution to the multispecies Lotka-Volterra equations, Revue SIAM (Society for Industrial and Applied Mathematics), 1994, Vol. 36(3) (1994), 480–488 https://doi.org/10.1137/1036104
    13. Ruan, J., Tan, Y., Zhang, C., A Modified Algorithm for Approximate Solutions of Lotka-Volterra systems, Procedia Engineering, 2011, Volume 15, 2011, Pages 1493-1497 https://doi.org/10.1016/j.proeng.2011.08.277
    14. Scarpello, G.M. and Ritelli, D., A new method for the explicit integration of Lotka--Volterra equations, 2003, Divulgaciones Matematicas, Vol. 11, No. 1, pp. 1--17.
    15. Seo, Gunog and Wolkowicz, Gail S. K., Sensitivity of the dynamics of the general Rosenzweig--MacArthur model to the mathematical form of the functional response: a bifurcation theory approach. Journal of Mathematical Biology, 76, No 7, 2018, 1873--1906.

    Return to Mathematica page
    Return to the main page (APMA0340)
    Return to the Part 1 Matrix Algebra
    Return to the Part 2 Linear Systems of Ordinary Differential Equations
    Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
    Return to the Part 4 Numerical Methods
    Return to the Part 5 Fourier Series
    Return to the Part 6 Partial Differential Equations
    Return to the Part 7 Special Functions


    J. D. Aplevich,Implicit Linear Systems, Lect. Notes in Control & Information Sci No. 152, Springer-Verlag, New York 1991.

    K. E. Brenan, S. L. Campbell and L. R. Petzold,Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, Elsevier, Amsterdam 1989.

    G. D. Byrne and W. E. Schiesser (eds),Recent Developments in Numerical Methods and Software for ODEs/DAEs/PDEs, World Scientific, Singapore 1991.

    S. L. Campbell,Singular Systems of Differential Equations, Pitman, London 1980.

    S. L. Campbell,Singular Systems of Differential Equations II, Pitman, London 1982.

    S. L. Campbell,The numerical solution of higher index linear time varying singular systems of differential equations, SIAM J. Sci. Stat. Comp. 334–348 (1985).

    S. L. Campbell,Rank deficient least squares and the numerical solution of linear singular implicit systems of differential equations, Contemp. Maths47, 51–63 (1985).

    S. L. Campbell,Consistent initial conditions for linear time varying singular systems, inFrequency Domain and State Space Methods for Linear Systems, Edited by C. I. Byrnes and A. Lindquist, Elsevier Sci Publ, Amsterdam 1986, pp. 313–318.

    S. L. Campbell,A general form for solvable linear time varying singular systems of differential equations, SIAM J. Math. Anal.18, 1101–1115 (1987).

    S. L. Campbell,A computational method for general higher index singular systems of differential equations, 1989 IMACS Transactions Sci. Computing1(2), 555–560 (1989).

    S. L. Campbell,Uniqueness of completions for linear time varying differential algebraic equations, Linear Algebra & Its Appl.161, 55–67 (1992).

    S. L. Campbell,Least squares completions for nonlinear differential algebraic equations, Numer. Math.65, 77–94 (1993).

    S. L. Campbell and C. W. Gear,The index of general nonlinear DAEs, preprint 1993.

    S. L. Campbell and E. Griepentrog,Solvability of general differential algebraic equations, SIAM J. Sci. Comp., to appear.

    S. L. Campbell and C. D. Meyer, Jr.,Generalized Inverses of Linear Transformations, Dover Press, New York 1991.

    S. L. Campbell and E. Moore,Constraint preserving integrators for general nonlinear higher index DAEs, Numer. Math., to appear.

    S. L. Campbell, E. Moore, and Y. Zhong,Utilization of automatic differentiation in control algorithms, IEEE Trans. Automatic Control39, 1047–1052 (1994).

    S. L. Campbell and W. J. Terrell,Observability of linear time varying descriptor systems, SIAM J. Matrix Analysis12, 484–496 (1991).

    S. L. Campbell, N. K. Nichols, and W. J. Terrell,Duality, observability, and controllability for linear time varying descriptor systems. Circuits Systems & Signal Process.10, 455–470 (1991).

    L. Dai,Singular Control Problems, Springer-Verlag, Berlin 1989.

    E. Griepentrog and R. März,Differential-Algebraic Equations and Their Numerical Treatment, Teubner-Texte zur Mathematik, Band 88, Leipzig 1986.

    E. Hairer, C. Lubich and M. Roche,The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods, Springer-Verlag, New York 1989.

    E. J. Haug and R. C. Deyo, Editors.Real-Time Integration Methods for Mechanical System Simulation, Springer-Verlag Computer & Systems Sci. Vol. 69, 1991.

    H. Krishnan and N. H. McClamroch,Tracking reference inputs in control systems described by a class of nonlinear differential algebraic equations, Proc. 1991 Conf. Dec. & Control.

    R. März,On quasilinear index 2 differential algebraic equations, Seminarberichte Nr. 92-1, Humboldt-Universität zu Berlin, Fachbereich Mathematik 1992, pp. 39–60.

    P. J. Rabier and W. C. Rheinboldt,A general existence and uniqueness theorem for implicit differential algebraic equations, Diff. Int. Eqns.4, 563–582 (1991).

    P. J. Rabier and W. C. Rheinboldt,A geometric treatment of implicit differential-algebraic equations, J. Diff. Eqns. (to appear).

    S. Reich,On a geometric interpretation of differential-algebraic equations, Circuits Systems & Signal Process.9, 367–382 (1990).

    S. Reich,On an existence and uniqueness theory for nonlinear differential-algebraic equations, Circuits Systems & Signal Process.10, 343–359 (1991).

    S. Reich,Symplectic integration of constrained Hamiltonian systems by Runge-Kutta methods, University of British Columbia Dept of Computer Sci. Techn. Rep. 93-13, 1993.

    S. Reich,On the local qualitative behavior of differential-algebraic equations, Circuits Systems & Signal Process, (to appear).

    C. Tischendorf,On stability of solutions of autonomous index-1 tractable and quasilinear index-2-tractable DAEs, Circuits Systems & Signal Process.13, 139–154 (1994).


    Math 216 Demonstrations

    This revisits the 5.2Whales demo, with explicit treatment of the linearization of the system. In scaled variables, we model the populations of krill ((x_1)) and baleen whales ((x_2)) with the system [ egin x_1' &= r_1 x_1 (1 - F_1 - x_1 - u x_2) x_2' &= r_2 x_2 left(1 - F_2 - frac ight), end ] where (r_1), (r_2), (F_1), (F_2) and ( u) are parameters. In this demonstration we take (F_1 = 0.5), (F_2 = 0), ( u = 1), in which case there are the equilibrium solutions ((x_1,x_2) = (1/4,1/4)) and ((x_1,x_2) = (1/2,0)). We consider (r_1 = r_2 = 2/5) and graph the solutions to the problem linearized at the equilibrium points along with trajectories found from the nonlinear system.

    Use Cases

    Lecture: The nonlinear equations may be presented with a minimal explanation of the different terms and used as an example of a system for which critical points and linear behavior there may be found. The demonstrations then graph these behaviors.

    Outside of Lecture: Solve the nonlinear system to find the critical points, and then find the linear systems approximating the nonlinear system at each. Show that the behavior that you see from the linear system is consistent with what the demonstrations show.

    Model Description

    As in 5.2Whales, we consider the interaction between a baleen whale species and its food-source, krill. Baleen whales feed by swimming through seawater in which krill (a crustacean about 5 centimeters long) is found with their mouths open, and then forcing the water out through its baleen, which filters the krill from the water.

    We consider the population of krill to be governed by its growth rate, constrained by environmental resources, and reduced by predation and fishing. Similarly, we take the population of whales to increase with a growth rate and be constrained by an environmental limiting factor inversely proportional to the krill population and reduced by predation.

    ODE Model

    Following 5.2Whales, the krill population (p_1) satisfies a modified logistic equation, [ p_1'(t) = r_1 p_1(t) left( 1 - frac ight) - C p_1(t) p_2(t) - r_1 F_1 p_1(t), ] where (p_2) is the whale population, so that the term (c p_1(t) p_2(t)) is the predation term (predation requires an interaction between the populations, and so is proportional to their product), and (r_1 F_1 p_1(t)) is the fishing term (the amount of krill caught is proportional to their population, and scaled as a fraction (F_1) of their growth rate).

    Similarly, if we assume that the carrying capacity for the whale population is inversely proportional to the krill population, as is suggested in [3], we obtain the equation [ p_2'(t) = r_2 p_2(t) left( 1 - frac ight) - r_2 F_2 p_2(t) ] for (p_2), where the term (r_2 F_2 p_2(t)) is again the fishing term. (This model is developed in [3].)

    We can nondimensionalize the populations (as in [3], [4], or [5]) by taking (x_1(t) = p_1(t)/K) and (x_2(t) = p_2(t)/(alpha K)), that is, by writing the populations as a fraction of their theoretical maxima. Introducing these variables and simplifying, we obtain the system [ egin x_1' &= r_1 x_1 (1 - F_1 - x_1 - u x_2) x_2' &= r_2 x_2 left(1 - F_2 - frac ight), end ] where ( u = Calpha K/r_1). Solving for the equilibrium solutions, we find there is a single non-zero equilibrium at [ x_1 = frac<1 - F_1><1 + u(1 - F_2)>,qquad x_2 = frac<(1 - F_1)(1 - F_2)><1 + u(1 - F_2)>. ] We will consider solutions near this equilibrium solution. There is a second equilibrium with (x_2 = 0), (x_1 = 1 - F_1 = 1/2), which we may consider as well.

    Near the equilibrium solution ((1/4,1/4)), and taking ( u = 1) (as in [3]), (r_2 = 0.4) (which is suggested by [6]), (F_2 = 0) (assuming that the whaling ban is actually observed) and the somewhat arbitrarily chosen value (F_1 = 0.5), we obtain the linearized system [ egin u_1' &= -frac14 r_1 u_1 - frac14 r_1 u_2 u_2' &= frac25 u_1 - frac25 u_2, end ] where (u_1) and (u_2) are the displacements from the equilibrium solution (in the scaled variables, ((1/4, 1/4))). We note that (u_1=0) and (u_2=0) corresponds to (x_1) and (x_2) being equal to their equlibrium values, (1/4). With (r_1 = 0.4) as well, the eigenvalues of the coefficient matrix are (lambda = -frac14pm ifrac><20>).

    Similarly, near ((1/2, 0)), we have the coefficient matrix (egin -r_1/2 & -r_1/2 0 & r_2end), and with (r_1 = r_2 = 0.4), eigenvalues are (lambda = -frac15, frac25).

    Matlab Demos

    Our demonstrations here show the solutions near each equilibrium solution, and those solutions in the full phase plane.

      Whales_Krill_Coexist.m: A demonstration that looks at the solutions near the coexistence point. The solutions to the linearized system are graphed in and the corresponding trajectory shown in the phase plane, along with a number of other representative trajectories. These are compared with the solution to the nonlinear problem, and then numerical solutions to the nonlinear system are shown for a number of initial conditions in the phase plane. Noter: also requires the file plot_localtraj.m. [show figure]

    Looking at the Model

    Some questions that may be worth considering:

    • What do the phase portraits near the equilibrium solutions tell us about the behavior of the system for the full phase plane?
    • How do we determine the direction of the spiral trajectories that occur when there are complex eigenvalues?

    References

    1. Wikipedia, Blue Whale. Wikipedia.org. Retrieved on: 23 Oct 2012
    2. Wikipedia, Fin Whale. Wikipedia.org. Retrieved on: 23 Oct 2012
    3. May, R.M., Beddington, J.R., Clark, C.W., Holt, S.J. and R.M. Laws (July 1979). "Management of Multispecies Fisheries." Science205(4403): 267-277.
    4. Edelstein-Keshet, L. Mathematical Models in Biology, SIAM Classics in Applied Mathematics 46. SIAM, 2005.
    5. Greenwell, R.N. Whales and Krill: A Mathematical Model, UMAP Module 610. COMAP, 1983.
    6. Beddington, J.R., and R.M. May (November 1982). "The Harvesting of Interacting Species in a Natural Ecosystem." Scientific American. November 1982: 62-69.