2. Méthodes numériques avancées et programmation#
Marc BUFFAT,dpt mécanique, Lyon 1
certains exemples sont tirés du cours « Scientific Python Lectures » de Robert Johansson
%matplotlib inline
from IPython.display import HTML,display
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 14
from matplotlib import animation
def printmd(string):
display(Markdown(string))
from metakernel import register_ipython_magics
register_ipython_magics()
2.1. Algorithmique et calcul scientifique sous Python#
2.1.1. Exigences du calcul scientifique#
« Le logiciel (software) est une des pierres angulaires de la science moderne. Sans logiciel, la science du vingt et unième siècle serait impossible. Mais sans de meilleurs logiciels, la science ne peut pas progresser »
importance de la documentation
importance de la validation
importance de la maîtrise du code (logiciels libres)
comprendre la méthode/démarche plus que la syntaxe
l’algorithme est plus important que le langage
# zen du Python
import this
import importlib
importlib.reload(this)
The Zen of Python, by Tim Peters
Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do it.
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!
The Zen of Python, by Tim Peters
Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do it.
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!
<module 'this' from '/usr/lib/python3.10/this.py'>
2.1.2. Démarche du calcul scientifique#
formulation mathématique du problème
recherche d’une solution (eventuellement approchée par discrétisation)
détermine ma méthode / algorithme pour le calcul de cette solution
top down design
programmation de cette méthode
bottom up programming (using libraries)
validation
analyse
2.2. Questions préliminaires#
2.3. Un premier exemple#
2.3.1. Formulation: calcul de la somme d’une série#
D.L. ordre N de \(log(1+x)\) $\( \lim_{N\rightarrow\infty} S_N = log(1+x) \)$
2.3.2. algorithmique#
1ere solution: traduction de la formule
Algorithme Serie(x,n)
import numpy as np
def Serie(x,n):
Sn = 0.
for i in range(1,n+1):
Sn = Sn + ((-1)**(i+1)) * (x**i) / i
return Sn
#
x0 = 0.2
S0 = Serie(x0,20)
print("S={} log(1+x)={}".format(S0,np.log(1+x0)))
S=0.18232155679395456 log(1+x)=0.1823215567939546
def Serie1(x,n):
Sn = 0.
coef = 1.
for i in range(1,n+1):
coef = -x*coef
Sn = Sn - coef / i
return Sn
x0 = 0.2
S1 = Serie1(x0,20)
print("S={} log(1+x)={}".format(S1,np.log(1+x0)))
S=0.18232155679395456 log(1+x)=0.1823215567939546
2.3.3. Solutions algorithmiques#
1ere solution: traduction de la formule
Algorithme Serie(x,n)
somme = 0
pour i de 1 a n
somme = somme + (-1)^(i+1)*x^i/i
retour somme
Solution non optimale !
2.3.4. Programmation Python#
écriture d’une fonction
validation
def Serie(x,n):
'''calcul serie S=x-x^2/2 + x^2/3 '''
somme = 0
return somme
# test
x0 = 0.2
S = Serie(x0,20)
# solution
def Serie1(x,n):
'''calcul serie S=x-x^2/2 + x^2/3 '''
somme = 0
for i in range(1,n+1):
somme = somme + ((-1)**(i+1)) * (x**i)/i
return somme
# test
x0 = 0.2
S1 = Serie1(x0,20)
print("S={} lim={}".format(S1,np.log(1+x0)))
S=0.18232155679395456 lim=0.1823215567939546
# version optimisée
def Serie2(x,n):
""" calcul de la somme de n termes de la serie x-x^2/2+x^3/3- """
coef=x
somme=0.0
for i in range(1,n+1):
somme = somme + coef/i
coef = -coef*x
return somme
# test
x0 = 0.2
S2 = Serie2(x0,20)
print("S={} lim={}".format(S2,np.log(1+x0)))
S=0.18232155679395456 lim=0.1823215567939546
%timeit Serie1(x0,1000)
369 μs ± 4.43 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
%timeit Serie2(x0,1000)
84.4 μs ± 497 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
sum([1,2,3])
[i**2 for i in range(20)]
[0,
1,
4,
9,
16,
25,
36,
49,
64,
81,
100,
121,
144,
169,
196,
225,
256,
289,
324,
361]
# version sans boucle
def Serie3(x,n):
return -sum([(-x0)**i/i for i in range(1,n+1)])
#
S3 = Serie3(x0,20)
print("S={} lim={}".format(S3,np.log(1+x0)))
S=0.18232155679395456 lim=0.1823215567939546
%timeit Serie3(x0,1000)
155 μs ± 636 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
2.4. Les langages de programmation scientifiques#
2.4.1. Langage machine#
suite d’instruction binaire executer par le processeur: 01010111 (suite 0 ou 1 : bits)
représentation en hexadécimal (octet ou bytes):
f8 = 248
instructions et données sont stockées en mémoire (à une adresse donnée)
@000141 instruction f8
spécifiques au processeur (intel, ARM, power, AMD, ..)
spécifique au système d’exploitation OS
%%bash
g++ serie.C
od -x -N 100 a.out
./a.out
0000000 457f 464c 0102 0001 0000 0000 0000 0000
0000020 0003 003e 0001 0000 1120 0000 0000 0000
0000040 0040 0000 0000 0000 3970 0000 0000 0000
0000060 0000 0000 0040 0038 000d 0040 001f 001e
0000100 0006 0000 0004 0000 0040 0000 0000 0000
0000120 0040 0000 0000 0000 0040 0000 0000 0000
0000140 02d8 0000
0000144
Calcul de la serie pour n=20 et x=0.2
Somme = 0.182322
Log(1+x)= 0.182322
2.4.2. Langage assembleur#
notation des instructions machines (move, push, load, add, mul)
programme + assembleur = code machine
%%bash
head -100 serie.s
.file "serie.C"
.local _ZStL8__ioinit
.comm _ZStL8__ioinit,1,1
.text
.globl _Z5seriedi
.type _Z5seriedi, @function
_Z5seriedi:
.LFB971:
.cfi_startproc
pushq %rbp
.cfi_def_cfa_offset 16
.cfi_offset 6, -16
movq %rsp, %rbp
.cfi_def_cfa_register 6
movsd %xmm0, -40(%rbp)
movl %edi, -44(%rbp)
movq -40(%rbp), %rax
movq %rax, -16(%rbp)
movl $0, %eax
movq %rax, -8(%rbp)
movl $1, -20(%rbp)
jmp .L2
.L3:
cvtsi2sd -20(%rbp), %xmm0
movsd -16(%rbp), %xmm1
divsd %xmm0, %xmm1
movapd %xmm1, %xmm0
movsd -8(%rbp), %xmm1
addsd %xmm1, %xmm0
movsd %xmm0, -8(%rbp)
movsd -16(%rbp), %xmm1
movsd .LC1(%rip), %xmm0
xorpd %xmm1, %xmm0
mulsd -40(%rbp), %xmm0
movsd %xmm0, -16(%rbp)
addl $1, -20(%rbp)
.L2:
movl -20(%rbp), %eax
cmpl -44(%rbp), %eax
jle .L3
movq -8(%rbp), %rax
movq %rax, -56(%rbp)
movsd -56(%rbp), %xmm0
popq %rbp
.cfi_def_cfa 7, 8
ret
.cfi_endproc
.LFE971:
.size _Z5seriedi, .-_Z5seriedi
.section .rodata
.LC3:
.string "Calcul de la serie pour n="
.LC4:
.string " et x="
.LC5:
.string "Somme de la serie = "
.LC7:
.string "Log(1+x)="
.text
.globl main
.type main, @function
main:
.LFB972:
.cfi_startproc
pushq %rbp
.cfi_def_cfa_offset 16
.cfi_offset 6, -16
movq %rsp, %rbp
.cfi_def_cfa_register 6
pushq %rbx
subq $56, %rsp
.cfi_offset 3, -24
movl %edi, -36(%rbp)
movq %rsi, -48(%rbp)
movl $20, -28(%rbp)
movabsq $4596373779694328218, %rax
movq %rax, -24(%rbp)
movl $.LC3, %esi
movl $_ZSt4cout, %edi
call _ZStlsISt11char_traitsIcEERSt13basic_ostreamIcT_ES5_PKc
movl -28(%rbp), %edx
movl %edx, %esi
movq %rax, %rdi
call _ZNSolsEi
movl $.LC4, %esi
movq %rax, %rdi
call _ZStlsISt11char_traitsIcEERSt13basic_ostreamIcT_ES5_PKc
movq %rax, %rdx
movq -24(%rbp), %rax
movq %rax, -56(%rbp)
movsd -56(%rbp), %xmm0
movq %rdx, %rdi
call _ZNSolsEd
movl $_ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_, %esi
movq %rax, %rdi
call _ZNSolsEPFRSoS_E
movl -28(%rbp), %edx
movq -24(%rbp), %rax
movl %edx, %edi
movq %rax, -56(%rbp)
2.4.3. Langage compilé#
C++ Fortran C
programme source + compilateur = programme binaire
programme source = portable
programme binaire = spécifique à l’ordinateur (non portable)
2.4.4. Langage C (Unix)#
%%bash
cat serie.C
#include <stdlib.h>
#include <math.h>
#include <iostream>
//
// calcul de la somme de n termes de la serie x-x^2/2+x^3/3-
//
double serie(double x,int n)
{
double coef=x;
double somme=0.0;
for (int i=1; i<=n; i++)
{
somme = somme + coef/i;
coef = -coef*x;
}
return somme;
}
//
int main(int argc,char *argv[])
{
int n=20;
double x=0.2;
std::cout<<"Calcul de la serie pour n="<<n<<" et x="<<x<<std::endl;
std::cout<<"Somme = "<<serie(x,n)<<std::endl;
std::cout<<"Log(1+x)= "<<log(1+x)<<std::endl;
return 0;
}
%%bash
g++ serie.C -o serie
./serie
Calcul de la serie pour n=20 et x=0.2
Somme = 0.182322
Log(1+x)= 0.182322
2.4.5. Langage interprété#
Python, Matlab
programme + interpréteur
portable = indépendant du système (Linux, MacOS, Windows)
Matlab
Matlab est un logiciel propriétaire qui est une boite à outils (toolbox) de calcul matriciel et numérique développé par Mathworks Inc. Il fonctionne sous Windows, Mac et Unix et permet d’écrire rapidement des scripts pour faire du calcul scientifique.
De nombreuses alternatives libres et de qualité existent, et notamment les logiciels: Scilab, Octave et surtout Python avec numpy.
%%bash
cat serie.m
% programme matlab
% ================
x=0.2; n=20;
% calcul de la serie
coef=x;
somme=0.0;
for i=1:n;
somme = somme+coef/i;
coef = -coef*x;
end;
display(sprintf('Calcul de la serie pour x=%g et n=%d',x,n))
display(sprintf('somme = %g',somme))
display(sprintf('log(1+x) = %g',log(1+x)))
exit
%%bash
# matlab -nojvm -r serie
Language Python
Python est un language de programmation moderne de haut niveau, logiciel libre, généraliste, multi-plateformes, multi-architecture, multi-OS (linux, windows, MacOS,..). Un programme python peut s’executer sur un tout petit raspberry Pi à 30€ (ARM), sur des smartphones, des portables, des PC jusqu’aux super-calculateurs HPC avec \(10^9\) coeurs de calcul.
 |
 |
 |
