L3 methodes numeriques 2015

Forum destiné aux étudiants de l'UGA (Université Grenoble-Alpes)

Modérateur : xcasadmin

Répondre
parisse
Messages : 5114
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

L3 methodes numeriques 2015

Message par parisse » lun. janv. 19, 2015 1:44 pm

Chapitre 1, representation des nombres sur machines
1/ Les entiers
- entiers courts 32 ou 64 bits, directement manipulables par le microprocesseur.
Rappel sur l'ecriture en base b par l'algorithme de division euclidienne repetee.
Limitations (calculs mod 2^32 ou 2^64)
- entiers longs: ecriture en base 2^32 ou 2^64 avec des entiers courts.
- rationnels.
Inconvenients: taille/temps de calcul -> autre type de representation

2/ Les nombres a virgule fixe: entier/denominateur fixe. Inconvenient: il faut beaucoup de chiffres pour pourvoir faire du calcul scientifique.
3/ Les nombres a virgule flottante (en abrege flottants): b fixe, n fixe +/- m*b^e avec m et e entiers, m ayant n chiffres en base b, e dans une plage [l,L] avec l<0 et L>0 et -l environ egal a L.
Flottant normalise: m>=b^(n-1), pour pouvoir avoir une representation la plus precise possible.
Erreur de representation: c'est une erreur relative<=b^(1-n) (et meme b^(1-n)/2 si on arrondit au plus proche)
Representation graphique "fractale" des nombres flottants.
Flottants denormalises: permet de faire le calcul si x!=y alors z=a/(x-y), lorsque x et y sont des flottants normalises dont l'exposant est (proche de) l'exposant minimal l.
Plus petit flottant (en-dessous on represente par 0), plus grand flottant, au-dela infinis. J'ai oublie de parler de NaN.
Les flottants sont des rationnels dont le denominateur divise une puissance de b, il est par exemple impossible de representer 1/10 en base 2, donc 0.5-5*0.1 ne donne pas forcement 0 (comme 2/3.-2*1/3. ne donne pas 0 en base 10).
Operations de base: par exemple mb^e*m'b^e'=mm'*b^(e+e'), mm' doit etre multiplie par b^-n ou b^(1-n) puis arrondi pour respecter le nombre de chiffres -> erreur relative d'arrondi.
L'erreur relative sur un produit est <= la somme des erreurs relatives + le produit (negligeable) + erreur d'arrondi finale. Donc environ 3*2^(-53) si 2^(-53) sur les operandes.
Exercice: pour l'inverse, au 1er ordre erreur relative sur l'argument+erreur de representation
L'erreur absolue sur une somme/difference est <= a la somme des erreurs absolues, ce qui peut faire perdre tous les chiffres significatifs du resultat si les arguments se compensent presque.
Calcul de f(x): erreur relative sur f(x) est sur |x*f'(x)/f(x)|, par exemple pour f(x)=exp(x), on multiplie par |x|
L'estimation de l'erreur sur par exemple la somme de N nombres avec N erreurs est pessimiste. Si on modelise les erreurs par des variables aleatoires d'esperance nulle et independantes, la variance de l'erreur totale est la somme des variances donc l'ecart-type est en sqrt(N) et pas en N.
Comment estimer l'erreur sur un calcul?
Estimation rigoureuse: difficile (a la main) et pessimiste. L'arithmetique d'intervalles permet d'automatiser ce type d'estimations.
Estimation empirique: on fait varier legerement les arguments et on regarde la sensibilite du resultat. Variante: on peut aussi utiliser des flottants multiprecision.
Ecriture scientifique des nombres flottants.

parisse
Messages : 5114
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: L3 methodes numeriques 2015

Message par parisse » lun. janv. 26, 2015 2:48 pm

26/1:
Retour sur NaN (not a number) oublie la derniere fois.

Types composes: complexes, listes/vecteurs/matrices/sequences.
Affectation d'un element de la liste,
* avec copie de liste := en temps O(n)
* affectation en place en temps O(1) mais a utiliser avec precaution!
Attention: si une variable t n'a pas ete initialise a une liste, t[1]:=3 cree une table et pas une liste.
Autres types plutot utiles pour faire du calcul formel:
* nom de variable, valeur affectee (ou hypothese assume()), purge
* expression
* fonction != expression

Exercice: calcul de la derivee numerique par les formules (f(x+eps)-f(x))/eps ou (f(x+eps)-f(x-eps))/eps, par exemple f(x):=sin(cos(tan(x))) et x=1.0. Comment choisir eps pour avoir une valeur la plus precise possible pour la derivee numerique?
Developpement de Taylor: f'(x)+O(eps) dans le 1er cas, f'(x)+O(eps^2) dans le second. Limitation: la precision relative est divisee par eps car les deux nombres que l'on retranche sont egaux a eps pres.
L'erreur est donc le max entre 2^-48/eps (precision par defaut en Xcas) et eps ou eps^2.
L'optimum est atteint lorsque les deux sont egaux, donc pour eps=2^-16 environ pour la 2eme formule et l'erreur est proche de 2^-32.

Chapitre 2 (debut) Algebre lineaire: autour du pivot de Gauss.
Rappel de l'algorithme: si A est une matrice LxC, on initialise l=1 et c=1,
tantque l<=L et c<=C
on cherche le pivot de valeur absolue/module maximal dans la colonne c a partir de la ligne l
si tous les coefficients sont nuls, on incremente c de 1 et on passe a l'iteration suivante du tantque
sinon on echange la ligne du pivot non nul avec la ligne l, puis on cree des 0 dans les lignes a partir de l+1 (cas (1) reduction sous-diagonale) ou pour toutes les lignes differentes de l (reduction complete).
on incremente l et c et on passe a l'iteration suivante du tantque.
Nombre d'operations (1 operation=1 multiplication et 1 addition) n^3/3 pour (1) et n^3/2 pour (2)
Applications
* resolution de systemes lineaires. On remarque que pour n grand, il vaut mieux faire la reduction sous-diagonale puis resoudre le systeme triangulaire que de faire la reduction complete.
* si AX=B et si on reduit la matrice A|B alors apres reduction on a A'X=B'
Par exemple si B=I et reduction complete A'X=B' avec A'=I donc X=B'=inverse de A. N.B.: algorithme pas tout-a-fait optimal
* recherche d'une base a partir d'une famille generatrice par reduction de la matrice avec les composantes de la famille generatrice en lignes
* recherche algorithmique d'une base reduite du noyau
* calcul du determinant

Interpretation matricielle de la reduction de Gauss: on suppose tout d'abord que le pivot est deja bien place (pas d'echange de lignes necessaires).
L'etape de creation de 0 colonne c revient a multiplier la matrice a reduire a gauche par une matrice L', obtenue a partir de l'identite en remplacant les 0 sous la diagonale en colonne c par les coefficients -alpha/p de la combinaison lineaire Lj <- Lj - alpha/p L_c (alpha coefficient ligne j colonne c, p pivot ligne c colonne c). L'inverse de L' s'obtient en changeant de signe ces coefficients (colonne c, lignes>c). Le produit des inverses de L' est une matrice L triangulaire inferieure, on a A=L*U, U etant la reduction sous-diagonale de A (matrice triangulaire superieure).
Je n'ai pas eu le temps de faire un exemple.

parisse
Messages : 5114
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: L3 methodes numeriques 2015

Message par parisse » mer. févr. 04, 2015 11:53 am

2/2:
Exemple 3,3 de factorisation LU sans echange de lignes pour [[1,2,3],[2,1,3],[2,0,1]], resolution d'un systeme lineaire en utilisant la factorisation LU.
Consequence des echanges de ligne, factorisation PA=LU
Applications: resolution de systemes lineaires, calcul d'inverse de matrice.
Interet en particulier si on doit resoudre plusieurs systemes avec meme matrice et differents seconds membres et que les seconds membres ne sont pas tous connus des le debut.

Cas des matrices symetriques definies positives: factorisation de Cholesky A=L*tran(L), calcul en posant le produit de matrices=A (cf. wikipedia) ou interpretation en termes de forme quadratique.
Exemple [[1,1,2],[1,6,1],[2,1,6]]
Generalisation: matrice hermitienne definie positive complexe A=LL*, factorisation A=LDL* avec L ayant des 1 sur la diagonale et D diagonale.

Rappel norme triple de matrice subordonnee a une norme sur L^n
Nombre de condition d'une matrice cond(A)=|||A|||*|||A^-1|||, lien avec l'erreur relative sur x en fonction de l'erreur relative sur b lorsqu'on resoud Ax=b
cond(A)>=1
Cas de la norme L2: cond(A)=plus grande valeur singuliere de A/plus petite valeur singuliere de A, ou les valeurs singulieres de A sont les racines carrees des valeurs propres (positives) de A*A.

Illustration avec une matrice de vandermonde
n:=5; H:=vandermonde(seq(j,j,1,n)); P:=egv(evalf(H)); w:=H^-1*tran(P)[0]; v:=H^-1*(tran(P)[0]+1e-3*tran(P)[4]); l2norm(w); l2norm(w-v)

Exercice: utiliser la decomposition lu de Xcas pour la matrice donnee en exemple au debut pour resoudre un systeme lineaire. Attention, il faut faire evalf() sur la matrice pour que la decomposition lu soit faite en prenant le pivot de module maximal sur une colonne.
Verifier que le temps de calcul pour lu est en O(n^3) sur une matrice aleatoire de taille n.
linsolve(p,l,u,b) sur les versions 1.3.22 a un bug il est en O(n^3) au lieu de O(n^2), bug corrige dans Xcas 1.4.3.

parisse
Messages : 5114
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: L3 methodes numeriques 2015

Message par parisse » lun. févr. 09, 2015 1:04 pm

9/2:
Interpolation polynomiale.
Polynome de Lagrange, demonstration de l'existence et unicite, ecriture dans la base de Newton 1, x-x0,..., (x-x0)*...*(x-x_{n-1}). Erreur f(x)-P(x) pour f fonction C^{[n+1]} en fonction de f^[n+1]. Majoration.
Choix des points d'interpolation equidistribues sur un intervalle [a,b], representation graphique de product(x-xi), plus grand aux bornes de l'intervalle. Phenomene de Runge: exemple pour 1/(1+x^2) sur [-10,10] avec n=20 et n=40. Augmenter le nombre de points n'est pas la solution, il vaudra mieux en choisir plus pres des bornes de l'intervalle.
Calcul efficace des coefficients dans la base de Newton par l'algorithme des differences divisees en O(n^2) operations.
Evaluation du polynome dans la base de Newton en 3n operations.
Sensibilite aux erreurs sur les f(x_i): l'erreur absolue est <= erreur absolue sur f(x_i)*somme(valeur absolue des polynomes de lagrange). On peut aussi voir l'interpolation comme la resolution d'un systeme lineaire de n+1 equations en n+1 inconnues (c'est inefficace car en O(n^3)) : plus le conditionnement de la matrice (de Vandermonde) est grand, plus les coefficients du polynome sont sensibles aux variations des f(x_i).

Exercices de programmation traites: differences divisees, et evaluation du polynome d'interpolation.

parisse
Messages : 5114
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: L3 methodes numeriques 2015

Message par parisse » mar. févr. 24, 2015 9:39 am

23/2: interpolation, points de Tchebyshev.
Polynomes de Tchebyshev T_n(cos(x))=cos(nx), relation de recurrence T_n+T_{n+2}=2XT_{n+1}, degre n, coeff dominant 2^(n-1) pour n>=1.
Propriete: |T_n(X)| <= 1 si |X|<=1, donc |produit_{i=0..n}(x-x_i)|<=1/2^n si les x_i sont les racines de T_{n+1} et si |x|<=1, majoration uniforme de la partie connue de l'erreur d'interpolation en les x_i.
Propriete: si U est de meme degre que T_n et de meme coefficient dominant, |U| ne peut etre strictement inferieur a 1 sur [-1,1]. Sinon T_n-U a le signe de T_n en les k*pi/n pour k=0..n, donc change de signe n+1 fois donc s'annule n fois, absurde car de degre<n.
Donc T_n est optimal pour la norme infinie sur [-1,1].
Si on choisit les x_i=(a+b)/2+(b-a)/2*t_i ou les t_i sont les racines de T_{n+1}, on a une erreur d'interpolation <= M_{n+1}/(n+1)!*1/2^n*((b-a)/2)^(n+1)
Interpolation par morceaux: pour approcher une fonction sur un grand intervalle [a,b] on utilise plusieurs polynomes, un par subdivision [alpha,beta], avec choix de alpha et beta assez proche pour que l'erreur ci-dessus soit faible.

Survol de quelques autres methodes d'approximation polynomiale:
- developpement de Taylor en un point x0, reste en (x-x0)^(n+1)/(n+1)!*f^{[n+1]}(xi_x)
Illustration de Taylor, Lagrange aux points de Tchebyshev, Lagrange en des points equidistribues pour sin(x) sur [-1,1]
- interpolation de Hermite, a mi-chemin entre Taylor et Lagrange, on impose la valeur de f et de derivees successives jusqu'a un certain ordre en certains points. Cela revient a faire tendre plusieurs points x_i vers une meme valeur dans l'interpolation de Lagrange. L'algorithme des differences divisees passe a la limite a condition de remplacer les quotients indetermines par la valeur des limites qui sont des derivees de la fonction au point.
- polynomes de Bernstein sur [0,1], B_n,k(x)=x^k*(1-x)^(n-k)*comb(n,k), proba de k succes si n tirages avec remise et x=proba de succes sur 1 tirage. Interet theorique pour montrer que les polynomes sont denses dans les fonctions continues sur [0,1] pour la norme infinie. Interet pratique pour n=3, courbes de Bezier avec 4 points de controle.
- approximant de Pade: fraction au lieu de polynome de Taylor, donne de meilleurs resultats, par exemple pour approcher la fonction exponentielle pres de 0.
- splines: approximation par intervalles avec conditions de recollement sur des derivees successives aux bornes de l'intervale

Exercices: approximation polynomiale de la fonction Si (sinus integral) sur [1,2] avec 11 points de Tchebyshev. Calcul de l'erreur theorique (en admettant une majoration de la derivee 11-ieme de Si sur [1,2]), l'erreur observee n'a plus aucun bit de mantisse significatif.

Interpolation de Hermite pour sin(x) sur [0,pi/2], en 0 et pi/2 on fixe la valeur de la fonction et de la derivee, calcul du polynome par l'algorithme des differences divisees, representation graphique de la fonction et du polynome sur [0,pi/2].

parisse
Messages : 5114
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: L3 methodes numeriques 2015

Message par parisse » lun. mars 02, 2015 2:32 pm

2/3: integration numerique
Subdivison, pas
Methodes classique: rectangles, trapezes, point milieu, Simpson
Exemple : calcul de area(x^3,x=0..1,10,..) selon les differentes methodes, observation de l'erreur.
Majoration de l'erreur pour rectangles M_1*h*(b-a)/2, trapezes M_2*h^2*(b-a)/12, point milieu M_2*h^2*(b-a)/24. Pour Simpson, on admet pour l'instant M_4*h^4*(b-a)/2880.
Nombre de subdivisions N pour avoir une erreur plus petite que epsilon: en 1/epsilon, 1/sqrt(epsilon) ou 1/epsilon^(1/4), environ 10^10 (long!), 10^5 (ca va), 10^2 (c'est mieux) pour epsilon=10^(-10). Les rectangles posent meme un probleme d'erreurs d'arrondis, a une precision de 10^(-15) il faudrait un jour pour calculer la valeur approchee, et l'erreur d'arrondi serait en 10^7!
Methode: I(f)=sum_j omega_j f(x_j) avec x_j points de la subdvision, omega_j de preference >0, sum omega_j=beta-alpha pour que la methode soit exacte.
Ordre d'une methode: entier n tel que la methode est exacte si appliquee a un polynome de degre <=n, mais ne l'est plus pour au moins un polynome de degre n+1.
Exemples: rectangle 0, trapezes et point milieu 1, Simpson au moins 2 (car interpolation en 2+1 points) en fait 3. On observe que la puissance de h dans l'erreur est ordre+1.
Une methode obtenue par interpolation en n+1 points est d'ordre au moins n.
Majoration universelle de l'erreur pour une methode d'ordre >=n en M_(n+1)*h^(n+1)/(n+1)!/2^(n+1)*(1+1/(n+2)) sous condition de positivite des omega_j (f de classe C_(n+1), M_(n+1) majorant de la derivee n+1-ieme en valeur absolue). Pour n=3, on trouve 1/320 comme constante (sub-optimale pour Simpson).
Majoration plus precise en utilisant le reste integral de la formule de Taylor a l'ordre n de f en alpha, noyau de Peano. Exemple cas de Simpson, on peut ainsi prouver la constante en 1/2880.
En pratique la valeur precise de la constante n'est pas tres importante sauf si on veut certifier une valeur approchee, on essaie plutot d'avoir une estimation de l'erreur pour utiliser une methode adaptative et faire varier la taille des subdivisions.
Optimisation de l'ordre par choix des points d'interpolations. Si on a n points d'interpolation, on peut avoir l'ordre 2n-1, en prenant les n racines du polynome de Legendre de degre n "rescale" de [-1,1] sur [alpha,beta], ou les polynomes de Legendre forment une famille orthonormale pour <f|g>=int(f(t)*g(t),t,-1,1)
Couple a une methode adaptative, on obtient ainsi une methode d'integration qui donne de tres bons resultats (cf. Hairer, estimation de l'erreur avec des quadratures emboitees).

parisse
Messages : 5114
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: L3 methodes numeriques 2015

Message par parisse » mar. mars 10, 2015 8:01 am

9/3: intégration (fin)
Accélération de la méthode des trapèzes: Richardson-Romberg.
Prop: si f est de classe C^2k, alors T_h(f) admet un développement limité à l'ordre 2k en puissances paires de h.
On peut éliminer la puissance h^2 en faisant T1_h(f)=(4*T_h/2(f)-T_h(f))/(4-1) (c'est la méthode de Simpson)
Puis la puissance h^4 en faisant (4^2*T1_h/2(f)-T1_h(f))/(4^2-1)
etc. on construit ainsi un tableau triangulaire, chaque ligne commence par les trapèzes sur une subdivision deux fois plus fine que la précédente, puis on applique les formules
(4^j*Tj-1_h/2(f)-Tj-1_h(f))/(4^j-1).
On s'arrếte lorsque la différence entre les deux derniers termes de 2 lignes consécutives est inférieure à la tolérance (erreur empirique).
Démonstration de la proposition: esquisse avec la formule d'Euler Mac Laurin et les polynômes de Bernoulli sur [0,1], puis sur [0,N] puis par changement de variable sur [a,b] découpé en N subdivisions de pas h=(b-a)/N.
Remarque: les coefficients des puissances paires de h du développement de T_h(f) font apparaitre en facteur la différence des dérivées impaires de la fonction en b et a, donc s'annulent si on intègre une fonction périodique sur une période, la méthode des trapèzes est donc très précise sans accélération (voir aussi l'exercice 5 de la feuille).

Autre estimation basée sur une méthode complètement différente: Monté-Carlo. Convergence en 1/sqrt(n) (admise) où n est le nombre d'évaluation. Illustration sur le calcul de int(exp(-t^2),t,0,2), estimation en prenant f(t):=exp(-t^2); n:=1000; 2*mean(apply(f,ranv(n,uniform,0,2)), histogramme sur plusieurs estimations avec m:=ranm(500,n,uniform,0,2); I:=seq(mean(f,m[k]),k,0,499); histogram(I,0,0.005);
l'écart-type est estimé empiriquement par stddevp(I) et en théorie par
sqrt(int(f(t)^2,t,0,2)/2-(int(f(t),t,0,2)/2)^2)/sqrt(n)

Exercices: programmation de Richardson-Romberg, calcul de l'ordre de l'exercice 4 (interpolation en les racines de Legendre pour n=2)

parisse
Messages : 5114
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: L3 methodes numeriques 2015

Message par parisse » lun. mars 16, 2015 1:26 pm

16/3: Methodes iteratives
1/ Point fixe
-> dans un intervalle ferme I de R:
definition de f contractante f:I-> I telle que |f(y)-f(x)| <= k*|y-x| avec k<1 pour x,y dans I
Thm: si f contractante alors u_{n+1}=f(u_n), u0 dans I, converge vers l'unique solution sur I de f(l)=l, et on a |u_{n+1}-l|<= |u_{n+1}-u_n|*k/(1-k), |u_n-l|<=C*k^n
Preuve faite avec les suites de Cauchy
Convergence lineaire (geometrique), pour avoir d digits le nombre d'iterations necessaire est proportionnel a d
Critere de contractance: recherche de majoration de |f'| sur I
Exemple f(x)=(x+2)/(x+1) sur I=[1,2], k=1/4, il faut a priori 17 iterations, en fait 14 suffisent pour avoir d=10 chiffres de l=sqrt(2), moitie moins pour 5 chiffres.
L'ecart vient du fait que quand on se rapproche de l, la constante de contractante peut etre majoree par un k proche de |f'(l)|
La vitesse de convergence n'est pas meilleure que geometrique si f'(l)!=0.

-> en dimension n (R^n, un ferme de R^n ou plus generalement un ferme d'un espace metrique complet)
L'enonce et sa preuve se generalisent.
Critere de contractance: plus difficile a etablir! On calcule la matrice de f' sur le ferme,
indice de ligne=numero de la composante de f
indice de colonne=variable de derivation
attention en Xcas, diff(f(x,y),[x,y]) donne la transposee de f'
Si la norme subordonnee de |f'| est <= k < 1 sur le ferme, alors f est contractante (on se ramene en dimension 1 en etudiant g(t)=f(x+t*(y-x))
Cas ou f est affine: f(x)=Cx+b, alors f'=C et on a juste une norme subordonnee a calculer.
Interet: resoudre plus rapidement des systemes de maniere approchee
on reecrit Ax=b en posant A=M-N, M diagonale ou en tout cas simple a inverser, sous la forme
x=f(x) avec f(x)=M^-1*(Nx+b)
Si k=||| M^-1*N ||| est <1, alors le point fixe converge, c'est particulierement interessant si k est petit.
Exemple: une matrice A 2*identite legerement perturbee, b dont toutes les composantes valent 1:
A:=2*idn(1000)+1e-4*ranm(1000,1000,uniform,-1,1) b:=[1$1000]; M:=diag(A)
Il faut 3 fois moins de temps pour trouver une solution approchee avec 10 decimales avec une methode de point fixe programmee a la main que par le pivot de Gauss. Pour une matrice 10000x10000, ce serait sans doute 30 fois moins!
Cette methode est la methode de Jacobi.
Critere de convergence: si A est a diagonale strictement dominante (|a_ii|>sum_{j != i} |a_{ij}|) alors la methode converge. Preuve: en calculant la norme subordonne a la norme infinie.
Autres methodes iteratives alternatives au pivot: Gauss-Seidel: on prend pour M la partie triangulaire inferieure, diagonale comprise, relaxation: M=1/omega*D+L (omega>0, D=diag(A), L=partie strictement sous la diagonale)
Critere de convergence: si A est symetrique definie positive et si transpose(M)+N l'est egalement, alors la methode converge. C'est en particulier le cas pour Gauss-Seidel si A est symetrique definie positive, et pour la relaxation si omega<2 et A symetrique definie positive.

Remarque: l'algorithme de classement des pages de google est une methode du point fixe (avec des donnees creuses, matrice de taille quelques milliards avec tres peu d'entrees non nulles).

2/ Methode de Newton (debut)
En dimension 1:
Interpretation geometrique: intersection de la tangente au graphe avec l'axe des abscisses
u_{n+1}=u_n-f(u_n)/f'(u_n)
Thm: si r est racine simple de f (f(r)=0, f'(r)!=0) et f de classe C^2 dans un voisinage [r-epsilon,r+epsilon] avec |f''|<=M2 et |1/f'|<=m1, il existe eta>0 tel que si |u_0-r|<eta alors la methode converge.
Preuve: faire Taylor a l'ordre 1 en u_n au point r pour obtenir |u_{n+1}-r|<=(M2*m1)/2*|u_n-r|^2
On a par recurrence la majoration |u_n-r|<= C^{2^n-1}*|u_0-r|^{2^n} avec C=(M2*m1)/2

Thm 2: si f est convexe sur I=[r,b] et f'(r)>0 alors la suite u_n de Newton converge en decroissant pour tout u_0 dans I.
(Demonstration au prochain cours).

parisse
Messages : 5114
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: L3 methodes numeriques 2015

Message par parisse » lun. mars 23, 2015 3:41 pm

23/3: evaluation de TP, un sujet tiré au sort parmi
* intégration numérique (N subdivisions, quadratures de Gauss à 7 points comparée avec interpolation en 7 points équidistribués)
* programmation de LU puis de la résolution des systèmes linéaires triangulaires inf/sup (on suppose pour commencer qu'il n'est pas nécessaire de permuter de lignes)
* programmation des méthodes de Jacobi, Gauss-Seidel et relaxation (exercice 10 de la feuille)
* recherche de toutes les racines approchées d'un polynôme à coefficients complexes (donné par la liste de ses coefficients).

parisse
Messages : 5114
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: L3 methodes numeriques 2015

Message par parisse » lun. mars 30, 2015 12:45 pm

30/3: methode de Newton (fin)
Demonstration de la convergence dans le cas convexe.
Estimation de l'erreur 0<=u_n-r<=f(u_n)/f'(r) par les accroissements finis.
Attention le calcul en flottants de f(u_n) fait generalement intervenir des nombres qui se compensent, f(u_n) peut renvoyer un resultat sans aucun chiffre correct en mantisse (voir negatif), il faut alors faire une evaluation en flottant multi-precision ou en calcul exact pour avoir une estimation plus precise de l'erreur.

Newton en dimension d
Theoreme: Soit f d'un ouvert U de R^d -> R^d de classe C^2, ayant une racine simple r (f(r)=0 et f'(r) inversible). Il existe une boule B(r,epsilon) telle que pour tout element u0 de B, la suite recurrente u_{n+1}=u_n-f'(u_n)^-1*f(u_n) converge vers r.
Ici f'(u_n) est la matrice de coefficient ligne i, colonne j df_i/dx_j.
Attention Xcas renvoie la transposee de f' si on calcule diff([f1,...,fd],[x1,...,xd])

Preuve analogue a la dimension 1, mais il faut utiliser le reste integral de Taylor.

Le calcul du rayon epsilon de la boule est en general impraticable. On peut utiliser Newton en dimension d soit en iterant avec l'espoir d'entrer dans un bassin d'attraction d'une racine, ou bien en complement d'une autre methode qui donne une valeur approchee grossiere de r, on affine alors la precision rapidement.

Exemple: algorithme de Durand-Kerner-Weierstrass pour calculer simultanement toutes les racines d'un polynome P: on developpe a l'ordre 1 en w Q(X,z+w)=product(X-(z_i+w_i))=P(X), on fait X=z_i ce qui donne le pas w_i=-P(z_i)/product_{j different de i}(z_i-z_j)
on part par exemple de z_j=exp(2*i*pi*j/d) ou d'une estimation des racines si on a une autre methode de recherche.
Proposition: il existe une racine de P dans la boule de centre z et de rayon degre(P)*|P(z)/P'(z)|

Exercice: implementation d'une fonction log si on dispose d'une fonction exp, en resolvant exp(z)=x en z

Répondre