3.6 Exercices

1.

Écrire une fonction pour résoudre un système linéaire où la matrice est triangulaire supérieure. On pourra utiliser l'instruction size qui permet de récupérer les deux dimensions d'une matrice :

-->[n,m]=size(A)

Dans un premier temps, on programmera l'algorithme classique utilisant deux boucles, puis on essaiera de remplacer la boucle interne par une instruction matricielle. Pour tester votre fonction, vous pourrez générer une matrice de nombres pseudo-aléatoires et n'en garder que la partie triangulaire supérieure avec l'instruction triu :

-->A=triu(rand(4,4))

2.

La solution du système d'équations différentielles du 1$^{\mbox{\scriptsize er}}$ ordre :

$$\frac{dx}{dt}(t) = A x(t), \quad x(0) = x_0 \in \mathbb{R}^n... ...(t) \in \mathbb{R} ^n, \quad A \in {\cal M} _{nn} (\mathbb{R})$$

peut être obtenue en utilisant l'exponentielle de matrice (cf votre cours d'analyse de 1$^{\mbox{\scriptsize ère}}$ année) :

x(t) = e^At x₀

Comme Scilab dispose d'une fonction qui calcule l'exponentielle de matrice (expm), il y a sans doute quelque chose à faire. On désire obtenir la solution pour $t \in [0,T]$. Pour cela, on peut la calculer en un nombre n suffisamment grand d'instants uniformément répartis dans cet intervalle $t_k = k \delta t, \quad \delta t = T/n$ et l'on peut utiliser les propriétés de l'exponentielle pour alléger les calculs :

$$x(t_k) = e^{A k \delta t} x_0 = e^{k (A \delta t)} x_0 = (e^{A \delta t})^k x_0 = e^{A \delta t} x(t_{k-1})$$

ainsi il suffit uniquement de calculer l'exponentielle de la matrice $A \delta t$ puis de faire n multiplications << matrice vecteur >> pour obtenir $x(t_1), x(t_2), \dots, x(t_n)$. Écrire un script pour résoudre l'équation différentielle (un oscillateur avec amortissement) :

$$x'' + \alpha x' + k x = 0, \mbox{ avec par exemple } \alpha = 0.1, k=1, x(0)=x'(0)=1$$

que l'on mettra évidemment sous la forme d'un système de deux équations du premier ordre. À la fin on pourra visualiser la variation de x en fonction du temps, puis la trajectoire dans le plan de phase. On peut passer d'une fenêtre graphique à une autre avec l'instruction xset("window",window-number). Par exemple :

--> //fin des calculs 
--> xset(``window'',0) // on selectionne la fenetre numero 0 
--> instruction pour le premier graphe (qui  s'affichera sur la fenetre 0) 
--> xset(``window'',1) // on selectionne la fenetre numero 1 
--> instruction pour le deuxieme graphe (qui s'affichera sur la fenetre 1)

3.

Écrire une fonction [i,info]=intervalle_de(t,x) pour déterminer l'intervalle i tel que $x_i \leq t \leq x_{i+1}$ par la méthode de la dichotomie (les composantes du vecteur x étant telles que x_i < x_i+1). Si $t \notin [x_1 , x_n]$, la variable booléenne info devra être égale à %f (et %t dans le cas inverse).

4.

Récrire la fonction myhorner pour quelle s'adapte au cas où l'argument t est une matrice (la fonction devant renvoyer une matrice p (de même taille que t) où chaque coefficient (i,j) correspond à l'évaluation du polynôme en t(i,j)).

5.

Écrire une fonction [y] = signal_fourier(t,T,cs) qui renvoie un début de série de Fourier en utilisant les fonctions :

$$f_1(t,T) = 1, \; f_2(t,T) = \sin(\frac{2 \pi t}{T} t), \; f... ...frac{4 \pi t}{T}), \; f_5(t,T) = \cos(\frac{4 \pi t}{T}),\cdots$$

au lieu des exponentielles. T est un paramètre (la période) et le signal sera caractérisé (en dehors de sa période) par le vecteur cs de ces composantes dans la base $f_1, f_2, f_3, \cdots$. On récupèrera le nombre de fonctions à utiliser à l'aide de l'instruction length appliquée sur cs. Il est recommandé d'utiliser une fonction auxiliaire [y]=f(t,T,k) pour calculer f_k(t,T). Enfin tout cela doit pouvoir s'appliquer sur un vecteur (ou une matrice) d'instants t, ce qui permettra de visualiser facilement un tel signal :

--> T = 1 // une periode ...  
--> t = linspace(0,T,101) // les instants ...  
--> cs =  [0.1 1 0.2 0 0  0.1] // un signal avec  une composante continue  
--> // du fondamental, pas d'harmonique  1 (periode 2T) mais une harmonique 2 
--> [y] = signal_fourier(t,T,cs); // calcul du signal 
--> plot(t,y) // et un dessin ...

6.

Voici une fonction pour calculer le produit vectoriel de deux vecteurs :

function [v]=prod_vect(v1,v2)
  // produit vectoriel v = v1 /\ v2
  v(1) = v1(2)*v2(3) - v1(3)*v2(2)
  v(2) = v1(3)*v2(1) - v1(1)*v2(3)
  v(3) = v1(1)*v2(2) - v1(2)*v2(1)

Vectoriser ce code de manière à calculer dans une même fonction function [v]=prod_vect_v(v1,v2) les produits vectoriels $v^i = v_1^i \wedge v_2^i$ où vⁱ, v₁ⁱ et v₂ⁱ désigne la i$^{\mbox{\scriptsize ème}}$ colonne des matrices (3,n) contenant ces vecteurs.