Informatique

Bonjour à tous,
C'est la première fois que je poste sur ce forum, et je suis ici pour demander de l'aide. Cela fait plusieurs mois déjà que je travaille à l'élaboration d'un code au format Python sur Jupyter Notebook. Je dois coder la simulation d'une voile solaire se déplaçant de la Terre vers Mars, et un second code dans lequel la voile se déplace de la Terre vers Vénus, en utilisant un RungeKutta d'ordre 4. Ainsi, en cliquant sur exécuter, je suis censé voir apparaître un graphique, sur lequel nous voyons le Soleil, la Terre et Mars représentés par des points de couleurs différentes, et un quatrième petit point ou autre chose quittant la Terre, pour aller vers l'orbite de Mars, là où cette dernière se trouvera pour l'intercepter.

Mars :
Vénus :

Pour ce faire, j'utilise les informations contenues dans ses deux documents :
http://www.nssc.ac.cn/wxzygx/weixin/201607/P020160718380095698873.pdf
https://www.researchgate.net/publication/353730296_Solar_sail_heliocentric_transfers_with_a_Q-law

Je commence à perdre espoirs, j'ai même au bout d'un moment cédé à ChatGPT pour qu'il m'aide à écrire un code, mais rien à faire. J'ai écris plusieurs codes, mais rien ne simule cela.
Voici mon code le plus abouti :

import numpy as np
import matplotlib.pyplot as plt

# Constantes
G = 6.6743e-11  # Constante gravitationnelle en m^3 kg^-1 s^-2
Ms = 1.989e30  # Masse du Soleil en kg
dMs = np.array([0, 0, 0])  # Position du Soleil
Mm = 6.39e23  # Masse de Mars en kg
dMm = np.array([2.067e11, 0, 0])  # Position de Mars
m = 10  # Masse de la voile solaire en kg
rho = 0.5  # Réflectivité de la voile solaire
A = 100  # Surface de la voile solaire en m^2
c = 299792458.0  # Vitesse de la lumière dans le vide en m/s
F0 = 1367.0  # Flux solaire à une distance moyenne de la Terre en W/m^2


# Fonction qui calcule la pression de radiation sur la voile solaire
def P(t, x, v):
    r = x - dMs
    F = F0 / (np.linalg.norm(r) ** 2)
    n = r / np.linalg.norm(r)
    theta = np.arccos(np.dot(-n, v) / np.linalg.norm(v))
    return (2 * F * np.cos(theta)) / c

# Fonction qui calcule l'accélération de la voile solaire
def a(t, x, v):
    Fv = P(t, x, v) * A * rho
    r = x - dMs
    Fg = -G * Ms * m * r / np.linalg.norm(r)**3
    Fm = -G * Mm * m * (x - dMm) / np.linalg.norm(x - dMm)**3
    return (Fv + Fg + Fm) / m

# Fonction qui calcule la trajectoire de la voile solaire avec la méthode de Runge-Kutta
def trajectoire(x0, v0, t0, tf, dt):
    t = [t0]
    x = [x0]
    v = [v0]
    while t[-1] < tf:
        k1x = v[-1]
        k1v = a(t[-1], x[-1], v[-1])
        k2x = v[-1] + k1v * dt/2
        k2v = a(t[-1] + dt/2, x[-1] + k1x * dt/2, v[-1] + k1v * dt/2)
        k3x = v[-1] + k2v * dt/2
        k3v = a(t[-1] + dt/2, x[-1] + k2x * dt/2, v[-1] + k2v * dt/2)
        k4x = v[-1] + k3v * dt
        k4v = a(t[-1] + dt, x[-1] + k3x * dt, v[-1] + k3v * dt)
        x_new = x[-1] + (k1x + 2*k2x +2*k3x + k4x) * dt/6
        v_new = v[-1] + (k1v + 2*k2v + 2*k3v + k4v) * dt/6
        t_new = t[-1] + dt
        x.append(x_new)
        v.append(v_new)
        t.append(t_new)
    return t, x, v


x0 = np.array([1.496e11, 0, 0]) # Position initiale de la voile solaire (Terre)
v0 = np.array([0, 30000, 0]) # Vitesse initiale de la voile solaire (Terre)
t, x, v = trajectoire(x0, v0, 0, 250000000, 100000)
x = np.array(x)
v = np.array(v)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x[:,0], x[:,1] ,'b')
ax.plot([dMs[0]], [dMs[1]], 'yo', markersize=10)
ax.plot([dMm[0]], [dMm[1]], 'ro', markersize=10)
ax.set_xlabel('Distance en m')
ax.set_ylabel('Distance en m')
ax.set_title('Trajectoire de la voile solaire vers Mars')
plt.show()

Il affiche ceci :

En résumé, je dois présenter un code ce Vendredi dans une UE (Mathématiques) coefficient 9 avec mon vieux code. Je sais que c'est un laps de temps un peu court, mais j'aimerais savoir si une âme charitable qui aurait déjà travaillé là-dessus ou non, serait prête à me partager un code simulant la trajectoire d'une voile solaire de la Terre vers Mars?
Je vous remercie de m'avoir lu,
Cordialement,

Et le probleme est?

Le problème c'est qu'en exécutant le code, il affiche quelque chose de statique qui ne bouge pas, c'est-à-dire cela au lieu de m'afficher quelque chose qui ressemblerait un peu plus à cela avec la ligne noire quittant la Terre qui se dessinerait progressivement jusqu'à atteindre le point représentant Mars.

Voici une bibliothèque de mécanique spatiale :

"""module keplerien bib_mclk, simplifiée en mk
period(a,mu)  retourne la periode en seconde
kepler(anomM,e):equation de kepler, donnée anomalie moyenne, resultat anomalie excentrique
keplerv(anomM,e):en radian equation de kepler, donnée anomalie moyenne, résultat anomalie vraie
invkepler(anomEx,e):equation de kepler inverse, donnée anomalie excentrique, résultat anomalie moyenne
invkeplerv(v,e): # en rd equation de kepler inverse,on donne l'anomalie vraie, résultat anomalie moyenne
kepCar(orb,mu):passe des kepleriens (a,e,i,pom,gom,v) aux cartésiens . dist en km, angles en rd  
carKep(pos,vit,mu):passage cartésien aux képlériens. Dist en km, angles en degré
extrapk (orb, ddatef,mu): #en rd ,date en seconde calcule l'anomalie vraie à la date +ddatef
gaussk(orbk,mu):matrice de gauss parametre képleriens
lamb(R1,R2,temps,mu)vecteur en majus, scalaire en minus, unité km et sec, methode lancaster limité aux ellipses et à moins d'un tour
"""

import math
import bib_math as mth
import numpy as np
import numpy.linalg as npl
"""ce module travaille sur une orbite képlérienne à 3 dimensions
il extrapole, calcule position, vitesse"""
           
mus = 1.327124e11  # km3/s2
                 
def period(a,mu):
    """retourne la periode en seconde"""
    per = 2*np.pi*math.sqrt(a*a*a/mu);
    return per
    
def kepler(anomM,e):
    """equation de kepler, donnée anomalie moyenne, resultat anomalie excentrique"""
    sinM = math.sin(anomM)
    cosM = math.cos(anomM)
    if (anomM ==0) or (anomM==np.pi):
        return anomM
    if (e==0):
        return anomM
    if (e<1):
        alpha = math.pow(math.pi,3)*e/np.pi;
        Z = e*sinM/math.sqrt(e*e+1-2*e*cosM);
        E0 = anomM+Z-Z*Z*Z*Z*cosM/(6*sinM*(1-e*e));
        if(E0<=0.7):
            E0 = anomM
        else:
            E0 = np.cbrt(anomM/alpha)
        for  n in range(6):
            sE0 = math.sin(E0)
            f = anomM-E0+e*sE0    
            fp = e*math.cos(E0)-1
            anomEx = E0-f/fp
            delta = f/fp   
            E0 = anomEx
      
    if (e>1):
        alpha = (-math.pi+math.pow(math.pi,3)*e)/math.pi
        E0 = math.cbrt(anomM*6/e)
        for n in range(6):
            sE0 = math.sinh(E0)
            f = anomM+E0-e*sE0    
            fp = -e*math.cosh(E0)+1
            fpp = -selfe*sE0
            anomEx = E0- f/fp
            delta = f/fp
            E0 = anomEx       
    return E0
         

def keplerv(anomM,e): # en radian
    """equation de kepler, donnée anomalie moyenne, résultat anomalie vraie"""
    anomEx = kepler(anomM,e);      
    if (e<=1):
        anomv = mth.angle(math.cos(anomEx)-e,math.sqrt(1-e*e)*math.sin(anomEx))
    if (e>=1):
        anomv = mth.angle(e-math.cosh(anomEx),math.sqrt(e*e-1)*math.sinh(anomEx))
    return anomv;

def invkepler(anomEx,e):
    """equation de kepler inverse, donnée anomalie excentrique, résultat anomalie moyenne"""
    if( e<=1):
        anomM = anomEx-e*math.sin(anomEx);
    if (e>=1):
        anomM =-anomEx + e*math.sinh(anomEx)
    return anomM;

     
def invkeplerv(v,e): # en rd
    """equation de kepler inverse, donnée anomalie vraie, résultat anomalie moyenne"""
    """angle en rd M fonction de v et ex"""       
    cosehf = math.cos(v)+e
    sinehf =math.sqrt(abs(1-e*e))*math.sin(v)  
    if (e<1):
        anomEx = mth.angle(cosehf,sinehf); 
        anomM = anomEx-e*math.sin(anomEx);
    if (e>1):
        anomEx = math.atanh(sinehf/cosehf);
        anomM = e*math.sin(anomEx)-anomEx;
    return anomM
     
def kepCar(orb,mu):
    """passe des kepleriens (a,e,i,pom,gom,v) aux cartésiens . dist en km, angles en rd """
    #indépendant de v   
    cpom = math.cos(orb[3])
    spom = math.sin(orb[3])
    cgom = math.cos(orb[4])
    sgom = math.sin(orb[4])      
    cosv = math.cos(orb[5])
    sinv = math.sin(orb[5])
    ci = math.cos(orb[2])
    si = math.sin(orb[2])
    P = [cpom*cgom-sgom*spom*ci,cpom*sgom+spom*cgom*ci,spom*si]
    Q = [-spom*cgom-cpom*sgom*ci,-spom*sgom+cpom*cgom*ci,cpom*si]

    p = orb[0] * (1 - orb[1] * orb[1]); #parametre p 
    #dépend de v
    r = p / (1 + orb[1] *cosv);
    racmusp = math.sqrt(mu/p);     
    F = r*cosv
    G = r*sinv; 
    pos = [0,0,0]
    for i in range(3):
        pos[i] = F*P[i] + G*Q[i]    
    FP = -racmusp*sinv
    GP = racmusp*(orb[1]+cosv);
    vit = [0,0,0]
    for i in range(3):
        vit[i] = FP*P[i] + GP*Q[i]
    return pos,vit

def carKep(pos,vit,mu):
    #print("type pos",type(pos))
    """passage cartésien aux képlériens. Dist en km, angles en degré"""
    pos = np.array(pos)
    #print("type pos",type(pos))
    vit = np.array(vit)
    orb = [0,0,0,0,0,0];
    r = npl.norm(pos,2)
    v = npl.norm(vit,2)
    #print("r,v",r,v)
    rv = np.dot(pos,vit)
    #print("rv ",rv)
    h = np.cross(pos,vit)
    hn = npl.norm(h)
   # print("h,hn ",h,hn)
    unsa = 2/r-v*v/mu;   
    if (unsa<=0):
        print("parabole")
        return
    p = hn*hn/mu;
    #print("p,mu ",p,mu)
    racpsmu = math.sqrt(p/mu);
    racmup = math.sqrt(mu*p);
    eCv = p/r-1;
    eSv = rv*racpsmu/r;
    #print("eCv,eSv",eCv,eSv)
    e2 = max(eCv*eCv+eSv*eSv,1e-18);
    #print("e2",e2)
    orb[5] = mth.angle(eCv,eSv);  #en rd
    U = [0,0,0]
    V = [0,0,0]
    for i in range(3):
        U[i] = pos[i]/r
        V[i] = (r*vit[i]-rv*U[i])/racmup
    #print("U, V ",U,V)
    l = mth.angle(U[0]+V[1], U[1]-V[0]);
    #print("l ",l*mth.crd)
    pompv = mth.crd*mth.angle(V[2], U[2]);  #en degré
    #print("pompv ",pompv*mth.crd)
    si2 = U[2]*U[2]+V[2]*V[2]
    orb[2] = math.asin(math.sqrt(si2))
    #print("si2,orb[2]",si2,orb[2])
    if orb[2]<=1e-10:
        orb[4] = 0
        pompv = l
    else:
        orb[4] = (l-pompv+2*np.pi)%(2*np.pi)  
        
    orb[1] = math.sqrt(e2)
    if orb[1]<1e-10:
        orb[3] = 0
        orb[5] = pompv
    else:
        orb[3] = (pompv-orb[5]+4*np.pi)%(2*np.pi)
    #orb[2] = pompv-orb[3] # en rd
    orb[0] = 1/unsa;
    return orb;

def extrapk (orb, ddatef,mu): #en rd ,date en seconde
    """calcule l'anomalie vraie à la date +ddatef en jours julien j2000. """
    result = orb[5];
    if (ddatef == 0):
        return result
    v=orb[5] #en rd"
    minitial = invkeplerv(v,orb[1]);   #M iniial en rd
    mn = math.sqrt(mu/math.pow(orbkep[0],3))# en rd/seconde    
    mfinal = (minitial+mn*ddatef)%(2*np.pi);# M final en rd   
    result = keplerv(mfinal,orb[1]);#v final tout en rd
   # print("\nv, minitial, mn, mfinal, result",v, minitial, mn, mfinal, result)  
    return result; # en rd
   
def gaussk(orbk,mu):
    """matrice de gauss parametre képleriens"""
    matgauss = np.zeros((6,3))
    sinv = math.sin(orbk[5])
    cosv = math.cos(orbk[5]);
    cosu = math.cos(orbk[5]+orbk[3])
    sinu = math.sin(orbk[5]+orbk[3]) 
    tgi = math.tan(orbk[2])      
    p = orbk[0]*(1-orbk[1]*orbk[1])
    racmup = math.sqrt(mu*p)
    racpsmu = math.sqrt(p/mu);
    r = p/(1+orbk[1]*cosv);
    matgauss[0][0] = 2*orbk[0]*orbk[0]*orbk[1]*sinv/racmup;   
    matgauss[0][1] = 2*orbk[0]*orbk[0]*racpsmu/r 
    matgauss[1][0] = sinv*racpsmu;
    aux = (orbk[1]+cosv)/(1+orbk[1]*cosv);
    matgauss[1][1] = racpsmu*(cosv+aux);
    matgauss[2][2] = r*cosu/racmup
    matgauss[3][0] = -racpsmu*cosv/orbk[1]; 
    matgauss[3][1] = racpsmu*sinv*(1+r/p)/orbk[1]; 
    matgauss[3][2] = -r*sinu/(racmup*tgi);
    matgauss[4][2] = r*sinu/(math.sin(orbk[2])*racmup)                                   
    matgauss[5][0] = - matgauss[3][0];
    matgauss[5][1] = - matgauss[3][1];
    matgauss[5][2] = 0;
    matgauss[6][0] = racmup/(r*r);
    return matgauss

def jul2000(date):
    """passe de la date civile aux jours julien ancien, wiki"""
    Y = +date[2];
    M = +date[1];  
    D = +date[0];  
    if (M<3):
       M = M+12
       Y = Y-1
    S = math.trunc(Y/100);
    B = 2-S+math.trunc(S/4);
    JJ = math.trunc(365.25*(Y+4716))+math.trunc(30.6001*(M+1))+D+B-1524-2451545;
    return JJ;

def cal2000(JJ):
    """passe des jours julien à la date civile jour, mois, année, heure//Wikipedia"""
    date = [0,0,0];
    JJ+=2451545;
    Z = math.trunc(JJ);
    F = JJ-Z;   
    a = math.trunc((Z-1867216.25)/36524.25);
    S = Z+1+a-math.trunc(a/4);  
    B = S+1524;
    C = math.trunc((B-122.1)/365.25);
    D = math.trunc(365.25*C);
    E = math.trunc((B-D)/30.6001);
    date[0] = B-D-math.trunc(30.6001*E)+F;
    if(E<13.5):
        date[1] = E-1
    else:
        date[1] =E-13 
    if(date[1]>2):
        date[2] = C-4716
    else:
        date[2] = C-4715
    return date;

def lamb(R1,R2,temps,mu):
    """vecteur en majus, scalaire en minus, unité km et sec"""
    """ methode lancaster limité aux ellipses et à moins d'un tour"""

    def lancaster(q,x):
        """x<1 ellipse"""
        K = q*q;"""K=1-c/s"""
        E = x*x-1;"""E<0, ellipse"""
        y = math.sqrt(-E);
        z = math.sqrt(1+K*E);
        l = mth.angle(x*z-q*E,y*(z-q*x));"""limité à un tour"""
        t = 2*(x-q*z-l/y)/E;
        dtsdx = (4-4*q*K*x/z-3*x*t)/E;
        return t,dtsdx,y,z
        
    code = 0;"""tout c'est bien passé"""
    R1=np.array(R1)
    R1é=np.array(R2)
    c = npl.norm(R2-R1);
    r1 = npl.norm(R1);"""aux.v;"""
    Ir1 = R1/r1;"""aux.U;"""
    r2 = npl.norm(R2);
    Ir2 = R2/r2; 
    #print("c,r1,Ir1,r2,Ir2",c,r1,Ir1,r2,Ir2)
    s = (r1+r2+c)*0.5;"""//demi-périmetre"""
    T=math.sqrt(8*mu/s)*temps/s; """//temps normalisé"""  
    IIh = np.cross(Ir1,Ir2);
    h = npl.norm(IIh);
    Ih = IIh/h;
    #print("s,T,IIh,h,Ih",s,T,IIh,h,Ih)
    teta = mth.angle(np.dot(Ir1,Ir2),h); """//cos(teta),sin(teta)"""
    #print("teta",teta);
    ctetas2 = M=math.cos(teta*0.5);
    q = math.sqrt(r1*r2)*ctetas2/s;
    if (Ih[2]<0):
        q =-q;
        aux = np.cross(Ir1,Ih);
        It1 = aux/npl.norm(aux);
        aux = np.cross(Ir2,Ih);
        It2 = aux/npl.norm(aux);
    else:
        aux = np.cross(Ih,Ir1) 
        It1 =aux/npl.norm(aux);
        aux = np.cross(Ih,Ir2)
        It2 = aux/npl.norm(aux);
        """approximation de la fonction par une hyperbole pour condition initiale"""
        T1 = 4*(1-q*q*q)/3;"""/pour x=1"""
        Tmin = T1*s/math.sqrt(8*mu/s);"""//temps en seconde pour x=1;"""
        """if (temps<Tmin):print("temps trop court. temps minimum = ", Tmin," temps=",temps);return};"""
        T0=2*(q*math.sqrt(1-q*q)+math.acos(q));"""//pour x=0, asymptote à x=-1"""
        xi= (T0-T)/(T+T0-2*T1); """ //condition initiale"""
       
        eps=1e-15;
        n=0;
        rl = lancaster(q,xi);
        #print("rl",rl)
        #print("rl(0)",rl[0])
        while ((n<10) and (abs(rl[0]-T)>eps)):  
            xi-=(rl[0]-T)/rl[1];
            rl = lancaster(q,xi);
            n +=1;
        """ //fin du while,"""
        if(n>20 and (abs(rl[0]-T)>eps)):
            print("pas de convergence,n=", n, "T-Ti =", T-rl.t); 
            code = 2;
            return code
        rl = lancaster(q,xi);
        #print("convergence n=",n," x=",xi," ti-T =",rl[0]-T);
        y = rl[2]; 
        z = rl[3];
        rp1=math.sqrt(2*mu*s)*(q*z*(s-r1)-xi*(s-r2))/(c*r1);
        rp2=math.sqrt(2*mu*s)*(-q*z*(s-r2)+xi*(s-r1))/(c*r2);
        unsa = 2*(1-xi*xi)/s;"""//2*y*y/s;"""
        e2 = (1-r1*unsa)*(1-r1*unsa)+r1*rp1*r1*rp1*unsa/mu;
        p = (1-e2)/unsa;      
        vteta1 = math.sqrt(mu*p)/r1;
        vteta2 = math.sqrt(mu*p)/r2;
        #print("vteta1,vteta2",vteta1,vteta2)
        V1 = Ir1*rp1 + It1*vteta1;
        V2 = Ir2*rp2 + It2*vteta2;
        #print("Lamb R1,V1,code,e2,p,unsa",R1,V1,code,e2,p,unsa);
        kep = carKep(R1,V1,mu);
        av2 = kep[5] + teta;
        return V1,V2,kep,av2,code
    """fin de lamb"""
    
    
if __name__ == "__main__": 
    
    print("\ntest de jul2000 et calend2000")
    jj = 40+10.5
    date = [31.65,10,2024]
    print("date",date)
    jj = jul2000(date)
    print("jj",jj)
    date = cal2000(jj)
    print("date",date)
    jj = jul2000(date)
    print("jj",jj)
    
    print("\ntest de lamb")
    ua = 150000000.0
    R1 = [ua,0,0]
    R2 = [0,ua,0]
    print("R1=",R1," R2=",R2)
    temps = 90*86400
    
    V1,V2,kep,av2,code = lamb(R1,R2,temps,mus)
    #print("\nV1,V2,kep,av2,code",V1,V2,kep,av2,code)
    orb1 = carKep(R1,V1,mus)
    print("\norb1",orb1)
    orb2 = carKep(R2,V2,mus)
    print("\norb2",orb2)
    
    print("\ntest de carKep/kepCar")
    R1 = [ua,0,0]
    V1 = [0,30.0,0]
    print("\nR1,V1",R1,V1)
    orbk = carKep(R1,V1,mus)
    print("\norbk",orbk)
    pos,vit = kepCar(orbk,mus)
    print("\npos,vit",pos,vit)
    pass 

Et voici une bibliothèque de Mathématique :

"""Module de mathématique.
modpositif(a,n) retourne le modulo positif
Vnorm(v) retourne la norme du vecteur à 3 dimension
angle(x,y) retourne l'angle polaire du vecteur (x,y) entre 0 et 2*pi
PScect(a,b) produit scalaire de a par b
PDTvect(a,b) produit vectoriel de a par b
quadGauss(npoints,a,b,f): retourne la quadrature de la fonction f entre a et b par la méthode de Gauss
RK4scal(t,y,h,f) intégre la fonction scalaire y'=f(t,y) sur le pas h par la méthode Runge-Kutta à l'orde 4 et retourne le nouvel y
RK4vect(t,y,h,f,nsys):intègr sur le pas h la fonction vectorielle y'=f(t,y) par la méthode Runge-Kutte ordre4 re retourne le nouvel y
zerosc(deb,fin,f,n):recherche de zero de la fonction f(x) entre a et b par dichotomie (section dorée à n coups)
zeroNewton(deb,f,dx,eps,nmax):zéro d'une fonction f(x) par la méthode de Newton avec dérivée numérique
    par jeux décalés de dx, jusqu'a la précision eps ou nmax itérations"""

import math

import numpy as np
import numpy.linalg as npl

pi = 3.141592653589793
pi = math.pi;
mut = 398601.3#km3/s2
mus = 132000000000.0#km3/s2
cdr = pi/180
crd = 180/pi

def modpositif(a,n):
    """donne le modulo de a par n toujours positif"""
    return ((a%n)+n)%n

def Vnorm(v):
    n = len(v);
    norm = 0;
    for i in range(n):
        norm += v[i]*v[i]
    return math.sqrt(norm)

def angle(lcos,lsin):
    """retourne l'angle entre 0 et 2*pi du vecteur x,y"""
    result = (math.atan2(lsin,lcos)+2*math.pi)%(2*math.pi)
    if result<0: 
        result+=2*pi 
    return result

def PSvect(a,b):
    """ utiliser plutôt np.dot"""
    return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]
    
def PDTvect(a,b):
    """utiliser plutôt np.coss"""
    c = [0,0,0]
    c[0] = a[1]*b[2]-a[2]*b[1];
    c[1] = -a[0]*b[2]+a[2]*b[0];
    c[2] = a[0]*b[1]-a[1]*b[0];
    return c

def quadGauss(npoints,a,b,f):
    """retourne la quadrature de la fonction f entre a et b par la méthode de Gauss"""
    A = (b-a)*0.5
    B = (b+a)*0.5
    pas=pi/npoints
    skk = [0 for i in range(npoints)]
    print(skk)
    quad = 0
    for i in range(npoints):#i va de zero à npoints-1
        sk = -math.cos((2*i+1)* pas*0.5)
        quad += math.sqrt(1-sk*sk)*f(A*sk+B)*pas
    quad *=A
    return quad
 
def RK4scal(t,y,h,f):
    """intégre la fonction scalaire y'=f(t,y) sur le pas h par la méthode Runge-Kutta à l'orde 4 et retourne le nouvel y"""
    k1 = h*f(t,y)
    k2 = h*f(t+h*0.5,y+k1*0.5)
    k3 = h*f(t+h*0.5,y+k2*0.5)
    k4 = h*f(t+h,y+k3)
    return y+(k1+2*k2+2*k3+k4)/6

def RK4vect(t,y,h,f,nsys):
    """intègr sur un pas h la fonction vectorielle y'=f(t,y) par la méthode Runge-Kutte ordre4 re retourne le nouvel y"""
    k1 = [0 for i in range(nsys)]
    k2 = [0 for i in range(nsys)]
    k3 = [0 for i in range(nsys)]
    k4 = [0 for i in range(nsys)]
    y1 = [0 for i in range(nsys)]
    y2 = [0 for i in range(nsys)]
    y3 = [0 for i in range(nsys)]
    yf = [0 for i in range(nsys)]
    k1 = f(t,y,nsys)
    for i in range(nsys):
        y1[i] = y[i]+h*k1[i]*0.5
    k2 = f(t+h*0.5,y1,nsys)
    for i in range(nsys):
        y2 [i]= y[i]+h*k2[i]*0.5
    k3 = f(t+h*0.5,y2,nsys)
    for i in range(nsys):
        y3[i] = y[i]+h*k3[i]
    k4 = f(t+h,y3,nsys)
    for i in range(nsys):
        yf[i] = y[i]+(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6      
    return yf
 
def zerosc(deb,fin,f,n):
    """recherche de zero de la fonction f(x) entre a et b par dichotomie """
    """"par section dorée à n coup"""
    nbdor = (1+math.sqrt(5))*0.5-1
    a = deb
    b = fin
    ya = f(a)
    if (ya ==0):
        return a    
    yb = f(b)
    if (yb == 0):
        return b
    if (ya*yb>0):
        return nan
    for i in range(n):
        c = a+(b-a)*nbdor
        yc = f(c)
        if ya*yc>0:
            a = c
            ya = yc
        else: 
            b = c
            yb = yc
    return (a+b)*0.5
    
def zeroNewton(deb,f,dx,eps,nmax):
    """zéro d'une fonction f(x) par la méthode de Newton avec dérivée numérique
    par jeux décalés de dx, jusqu'a la précision eps ou nmax itérations"""
    a = deb
    n = 0
    ya = f(a)
    while (abs(ya)>eps) and (n<nmax):
    #for i in range(10):
        yb = f(a+dx)
        dya = (yb-ya)/dx
        a -= ya/dya
        ya = f(a)
        n += 1
    return a
        



 
if __name__ == "__main__":
    
  
    print("pi=",pi)
    print("mut=",mut)
    print("mus=",mus)
    print("cdr=",cdr)
    print("cdr=",crd)
    
    print ("test de modpositif")
    a=-7
    b=5
    print(a,b,modpositif(a,b))
    
    print ("\n test de Vnorm")
    v = [2,3,4]
    norm1 = Vnorm(v)
    print("v, norme de v ",v,norm1)
    v = np.array([2,3,4])
    print("v, norme de v avec np",npl.norm(v,2))
   
    print("\ntest de ps")
    a = [1,3,5];
    b = [6,8,9]
    print("ps de a,b",PSvect(a,b))
    print("avec np.dot",np.dot(a,b))
    
    print("\ntest de PDTvect")
    a = [1,3,5];
    b = [6,8,9]
    print("PDTvect de a,b",PDTvect(a,b))
    print("avec np.cross",np.cross(a,b)) 
    
    print("\ntest de multpilcation par scalaire")
    a = [0,1,2,3,4]
    l = 2
    print("a,l et a*l",a,l,np.multiply(a,l))
    
    print('\n test de angle')
    a = math.cos(30*cdr)
    b = math.sin(30*cdr)
    c = angle(a,b)*crd
    print(a,b,c)
    
    print ("\n","test de quadGauss")
    def f(var):
        return var+1
    a=0
    b=1
    np=10
    c = quadGauss(np,a,b,f)
    print("quad=","%3.2f" % c,"qvrai=",0.5*b*b+b)
    
    print("\n","test de RK4scal")
    t=0
    y=0
    h=1
    def ff(t,y):
        return t*t+3*t
    c = RK4scal(t,y,h,ff)
    print("c=","%3.2f" % c)
    
    print("\n","test de RK4vect")
    t=0 
    h=1
    nsys=2
    y=[0 for i in range(nsys)]
    def fff(t,y,nsys):
        yp = [0 for i in range(nsys)]
        for i in range(nsys):
            yp[i] = t*t
        return yp
    c =  RK4vect(t,y,h,fff,nsys) 
    print("c=", c)
    
    print("\n test de zerosc(deb,fin,f,n)")    
    deb = -1
    fin = 3
    n = 10
    def fz(x):
        return x-1
    c = zerosc(deb,fin,fz,n)
    print("c ","%3.2f" % c)
    
    print("\n test de la fonction zeroNewton(deb,f,dx,eps,nmax)")
    deb = -1
    fin = 3
    dx = 1e-3
    nmax = 10
    eps = 1e-5
    def fz(x):
        return x-1
    c = zeroNewton(deb,fz,dx,eps,nmax)
    print("c= ","%3.2f" % c)

Dans l'attente que quelqu'un me réponde favorablement.

Le 18 mai 2023 à 17:53:57 :

Ça ne sert à rien de up si t'es déjà le premier topic de la liste.