Bonjour,
je suis actuellement sur un projet qui consiste à simuler et tracer en 3d la trajectoire d'un objet(vaisseau,comète,...) dans le système solaire. J'avais déjà fait le script lorsque l'objet n'était attiré que par un seul astre et ça marche très bien, aucun problème.
Mais je veux simuler le mouvement de l'objet dans le système solaire il faut donc que je rajoute toutes les autres planètes. J'ai donc rajouté la Terre(initialement il n'y avait que le soleil placé sur l'origine du repère). Et j'ai donc refait les équations mais cette fois avec 2 astres, et on peut donc voir qu'il y a un 2ème membre qui apparaît dans l'équation différentielle (j'ai fait le calcul avec les coordonnées cartésiennes) :

avec ici r=r1:

et:

Et donc les M sont les masses des astres, les r les distances entre le centre de l'astre i et l'objet.
x,y,z sont les coordonnées de l'objet dans le repère dont le centre du soleil est l'origine. Et x2,y2,z2 sont les coordonnées de la terre dans ce même repère.
Et N est le nombre d'astre donc ici N=2.
Et donc s'il n'y a qu'un seul astre, B est la matrice nulle.
Donc lorsqu'il n'y avait qu'un astre, B n'était pas présent et il n'y avait aucun problème.

clear
clc
 
G=6.67*10^-11;
//MASSES (kg)
M_Terre=5.98*10^24;//Masse de la Terre (kg)
M_Soleil=1.989*10^30;
//RAYONS (m)
R_Terre=6.371*10^6;//Rayon de la Terre (en m)
R_Soleil=696340*10^3;
//COORDONNEES DE CHAQUE ASTRE (sous la forme (x,y,z))
//Soleil c'est (0,0,0)
coordTerre=[149597870*10^3 0 0];//Distance terre-soleil puisque l'origine du repère est le soleil
 
function u_point=f(t,u)//vecteur u_point régissant les equations du mouvement
    //Norme des vecteurs r avec pour origine chaque astre
    r=sqrt(u(1)^2+u(2)^2+u(3)^2);//Norme r ayant pour origine l'origine du repère donc le premier astre(soleil)
    r2=sqrt((u(1)-coordTerre(1))^2+(u(2)-coordTerre(2))^2+(u(3)-coordTerre(3))^2);//r2 en partant de la terre
 
    M=[M_Soleil M_Terre];//Vecteur des masses
    r_cube=[r^3 r2^3];//Les rayons puissance 3
 
    vect_somme=M./r_cube;
 
    K=-G*sum(vect_somme);
 
    if r<R_Soleil then
        u_point=[0;0;0;0;0;0];
    elseif r2<R_Terre then
        u_point=[0;0;0;0;0;0];
    else
        A=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;K 0 0 0 0 0;0 K 0 0 0 0;0 0 K 0 0 0];//matrice A de l'équation u_point=A*u+B
        B=[0;0;0;G*M_Terre*coordTerre(1)/r2^3;G*M_Terre*coordTerre(2)/r2^3;G*M_Terre*coordTerre(3)/r2^3];//Ici c'est une somme mais elle n'est pas écrite car il n'y a que le soleil et la terre, il faudra faire la somme si il y a un 3eme astre
 
        u_point=A*u+B;//Utiliser une autre fonction à cause de B
    end
endfunction
 
function [x,y,z]=sphere(phi,teta)//calcul de la sphere
    x=cos(phi).*cos(teta);
    y=cos(phi).*sin(teta);
    z=sin(phi);
endfunction
 
function f=afficherAstre(origine,rayon)
    //origine est un vecteur [x y z] representant les coordonnées de l'astre
    //rayon est une valeur donnant le rayon de l'astre
 
    [x,y,z]=eval3dp(sphere,linspace(-%pi/2,%pi/2,40),linspace(0,2*%pi,20));//calcul des facettes de la sphere
    //on va positionner et dimensionner notre astre
    [nbLignes_x,nbColonnes_x]=size(x);
    x=x+origine(1)*ones(nbLignes_x,nbColonnes_x);//on positionne l'origine de la sphere
    [nbLignes_y,nbColonnes_y]=size(y);
    y=y+origine(1)*ones(nbLignes_y,nbColonnes_y);//on positionne l'origine de la sphere
    [nbLignes_z,nbColonnes_z]=size(z);
    z=z+origine(1)*ones(nbLignes_z,nbColonnes_z);//on positionne l'origine de la sphere
    //on dimensionne notre astre pour avoir le bon rayon
    x=x.*rayon;
    y=y.*rayon;
    z=z.*rayon;
 
    clf();
    plot3d(x,y,z);
endfunction
 
function u=rotation(distance,v0,hours)
    //distance est la distance entre l'objet et l'astre (en km)
    //v0 est la vitesse initiale (m/s)
    //hours est le temps de la simulation en heures
    u0=[R_Soleil+distance; 0; 0;v0/sqrt(3);v0/sqrt(3);v0/sqrt(3)];//vecteur v0 initial de la forme [x y z x_point y_point z_point]
    distance=distance*1000;
    t=linspace(0,3600*hours,2000);
    u=ode(u0,0,t,f);//resolution du systeme differentiel
 
    afficherAstre([0 0 0],R_Soleil);//affichage de l'astre
 
    //tracer la trajectoire
    comet3d(u(1,:),u(2,:),u(3,:),"colors",3);
endfunction
 
//DEBUT DU PROGRAMME
 
//Parametres initiaux
altitude=10^10;//en km sans compter le rayon du soleil donc c'est bien l'altitude
vitesse=20000;
duree=1000;
rotation(altitude,vitesse,duree);

Et voilà l'erreur:

Pour trouver une solution j'ai modifié mon code en ajoutant printf('%d %d - %d %d - %d %d\n', size(u), size(A), size(B)) juste avant la ligne u_point=A*u+B;//Utiliser une autre fonction à cause de B mais il se passe quelque chose de bizarre :

Le nombre de colonnes de u change d'un coup.

Si vous avez une idée je suis preneur.

Tu n'utilises pas le vecteur temps (t) dans ta fonction f, il doit y avoir quelque chose de "mal fait" avec lui car c'est le seul qui a une taille de 2000 points

Non justement dans les équations du mouvements t n'apparaît pas. De plus ça marchait très bien lorsqu'il n'y avait que A peut être que le problème viens de B mais je ne vois pas où

étonnant que le temps n'intervienne pas quand même

Enfaites pas tellement puisque la force c'est -G*m*M/r²

Sinon, je viens de tester le code sur https://cloud.scilab.in/ (après avoir supprimer les lignes vides et refait la tabulation). Je n'ai pas eu d'erreurs. et le résultat obtenu me donne l'image suivante:

Ah ouais, c'est bizarre. Je sais pas d'où ça vient alors.