Problème "horzcat" en utilisant ODE45

Fermé
mattard54 Messages postés 2 Date d'inscription samedi 1 juin 2013 Statut Membre Dernière intervention 14 juin 2013 - Modifié par mattard54 le 14/06/2013 à 11:07
Bonjour,

je débute sur matlab et je dois, dans le cadre de mon stage, écrire un programme qui me permettra de visualiser les trajectoires de particules chargées à l'intérieur d'un réacteur. Le problème vient du fait que matlab m'écrit
"Error using horzcat
CAT arguments dimensions are not consistent."
quand je lance mon programme. SI vous pouviez m'indiquer la raison pour laquelle ça ne marche pas ça m'aiderait vraiment.

Voici le programme que j'ai écris:

clear all;
global coeff_q R0 B0 m a q N X Y Z Exc Eyc Ezc theta phi rp Ert Ept Ett Ettot

% Petit rayon du tokamak
a=2.0;

% Grand rayon du tokamak
R0=6.2;

% Température
T=1.0e4*11600;

% Temps total d'intégration
ttot=4.0e-4;

% Coeff_q pour le profil de sécurité q(r) 
coeff_q=2.0;

% Nombre de particules dans le tokamak
N=2;

% Masse des particules
m= [9.1e-31, 3.3437e-27];

% Charge des particules
q= [-1.6e-19, 1.6e-19]; 

% Champ magnétique B_phi0 (en r=0)
B0=5.0;

% Rayon de départ des centres-guides
r_dep= [0.8*a, 0.79*a];

% Champs magnétiques toroïdals là où démarre les centres-guides
Bdep= zeros(N,1);
for i=1:N
Bdep(i)=B0*R0/(R0+r_dep(i));
end

% facteurs de sécurité là où démarre les centres-guides
q_secur= zeros(N,1);
for i=1:N
q_secur(i)=1.0+coeff_q*r_dep(i)^2/a^2;
end

% Norme des champs magnétiques là où démarre les centres-guides
B= zeros(N,1);
for i=1:N
B(i)=sqrt(Bdep(i)^2+(r_dep(i)*Bdep(i)/(R0+r_dep(i))/q_secur(i))^2);
end

% Pulsations cyclotron et périodes cyclotron là où démarre les
% centres-guides
periode= zeros(N,1);
omega_c= zeros(N,1);
for i=1:N
omega_c(i)=q(i)*B(i)/m(i);
periode(i)=2.0*pi./omega_c(i);
end

% Vitesse thermique des particules
vT= zeros(N,1);
for i=1:N
vT(i)=sqrt(1.38e-23*T/m(i));
end

% Rayon de Larmor là où démarre les centres-guides
rl= zeros(N,1);
for i=1:N
rl(i)=vT(i)/omega_c(i);
end

% rapport v_phi0 / v_S0 (presque v//0 / v_perp0) 
rapport=0.3;
v_S0= zeros(N,1);
v_phi0= zeros(N,1);
for i=1:N
v_S0(i)=sqrt(1.0-rapport)*vT(i);
v_phi0(i)=sqrt(rapport)*vT(i);
end

%Coordonnées cartésiennes des particules
X= [0.1, 0.2];
Y= [6, 6.1];
Z= [0.1, 0.2];

%Détermination de la distance entre chaque particules en cartésien
r= [i, j];
for i=1:N
    for j=1:N
        r(i,j)=sqrt((X(j)-X(i))^2+(Y(j)-Y(i))^2+(Z(j)-Z(i))^2);
    end
end

%Détermination du vecteur Rij en cartésien
rx= [i, j];
ry= [i, j];
rz= [i, j];
for i=1:N
    for j=1:N
        rx(i,j)=X(j)-X(i);
        ry(i,j)=Y(j)-Y(i);
        rz(i,j)=Z(j)-Z(i);
    end
end

%Calcul du vecteur directeur entre chaque particule en cartésien
ux= [i, j];
uy= [i, j];
uz= [i, j];
for i=1:N
    for j=1:N
        ux(i,j)=rx(i,j)/r(i,j);
        uy(i,j)=ry(i,j)/r(i,j);
        uz(i,j)=rz(i,j)/r(i,j);
    end
end

%Déterminatation de E pour chaque particule en cartésien
Exc= zeros(N,1);
Eyc= zeros(N,1);
Ezc= zeros(N,1);
Ex= [i, j];
Ey= [i, j];
Ez= [i, j];
for i=1:N
    Exc(i)=0;
    Eyc(i)=0;
    Ezc(i)=0;
    for j=1:N
       Ex(i,j)=(1/(4*pi*8.85e-12))*(q(j)/(r(i,j)^2))*ux(i,j);
       Ey(i,j)=(1/(4*pi*8.85e-12))*(q(j)/(r(i,j)^2))*uy(i,j);
       Ez(i,j)=(1/(4*pi*8.85e-12))*(q(j)/(r(i,j)^2))*uz(i,j);
       
       Exc(i)=Exc(i)+Ex(i,j);
       Eyc(i)=Eyc(i)+Ey(i,j);
       Ezc(i)=Ezc(i)+Ez(i,j);
    end
end


%Coordonnées toroïdales des particules

theta=[90, 89];
phi=[45, 46];
rp=[a/2, a/2.1];


%Ouverture de la boucle temporelle pour la résolution des équations
%différentielles à chaque pas de temps
t= zeros(201,1); %je prends 200 pas de temps pour arriver jusque ttot
t(1)=0;
for d=1:201
t(d+1)=t(d)+2e-6;
    
%Détermination de X Y et Z en coordonnées toroïdales
Xt= [i, j];
Yt= [i, j];
Zt= [i, j];
for i=1:N
    Xt(i)=(R0+rp(i)*cos(theta(i)))*cos(phi(i));
    Yt(i)=(R0+rp(i)*cos(theta(i)))*sin(phi(i));
    Zt(i)=-rp(i)*sin(theta(i));
end

%Détermination du vecteur Rij en toroïdale
rtr= [i, j];
rtt= [i, j];
rtp= [i, j];
for i=1:N
    for j=1:N
        rtr(i,j)=(X(j)-X(i))*(cos(theta(i))*cos(phi(i)))+(Y(j)-Y(i))*((Z(i)*cos(theta(i))/rp(i)^2))-(Z(j)-Z(i))*(Y(i)/(X(i)^2+Y(i)^2));
        rtt(i,j)=(X(j)-X(i))*(cos(phi(i))*sin(theta(i)))+(Y(j)-Y(i))*((Z(i)*sin(phi(i))/(rp(i)^2)))+(Z(j)-Z(i))*(X(i)/(X(i)^2+Y(i)^2));
        rtp(i,j)=(X(j)-X(i))*sin(phi(i))+(Y(j)-Y(i))*cos(theta(i))/rp(i);
    end
end

%Détermination de la distance entre chaque particule en toroïdale
rt= [i, j];
for i=1:N
    for j=1:N
        rt(i,j)=sqrt((Xt(j)-Xt(i))^2+(Yt(j)-Yt(i))^2+(Zt(j)-Zt(i))^2);
    end
end

%Détermination du vecteur directeur entre chaque particule en toroïdale
ur= [i, j];
ut= [i, j];
up= [i, j];
for i=1:N
    for j=1:N
        ur(i,j)=rtr(i,j)/rt(i,j);
        ut(i,j)=rtt(i,j)/rt(i,j);
        up(i,j)=rtp(i,j)/rt(i,j);
    end
end

%Détermination de E pour chaque particule en Toroïdale
Ert= zeros(N,1);
Ett= zeros(N,1);
Ept= zeros(N,1);
Ettot= zeros(N,1);
Er= [i, j];
Et= [i, j];
Ep= [i, j];
for i=1:N
    Ert(i)=0;
    Ett(i)=0;
    Ept(i)=0;
    Ettot(i)=0;
    for j=1:N
        Er(i,j)=((1/(4*pi*8.85e-12))*(q(j)/(rt(i,j)^2))*ur(i,j));
        Et(i,j)=((1/(4*pi*8.85e-12))*(q(j)/(rt(i,j)^2))*ut(i,j));
        Ep(i,j)=((1/(4*pi*8.85e-12))*(q(j)/(rt(i,j)^2))*up(i,j));
        
        Ert(i)=Ert(i)+Er(i,j);
        Ett(i)=Ett(i)+Et(i,j);
        Ept(i)=Ept(i)+Ep(i,j);
        Ettot(i)=Ert(i)+Ett(i)+Ept(i);
    end
end

% Résolution avec ode45 des 4 équations différentielles du 1er ordre
options = odeset('RelTol',1e-8,'AbsTol',[1e-16 1e-16 1e-16 1e-16]);
[T,X] = ode45(@trajCHANGE,[t(d) t(d+1)],[r_dep 0.0 0.0 v_phi0],options);

end





Et voici ce que j'ai écris dans la fonction trajCHANGE:

function dx = trajCHANGE(t,x);

% Fonction qui retourne le vecteur x en fonction du temps t 

global  coeff_q R0 B0 m a q N Btheta Bphi theta rt

% Rq : B0 est B_phi_0

% Calcul de R et B en fonction de la position du centre-guide
B= zeros(N,1);
R= zeros(N,1);
Btheta= zeros(N,1);
Bphi= zeros(N,1);
for i=1:N
R(i)=(R0+rt(i)*cos(theta(i)));
B(i)=sqrt(B0^2*R0^2/R(i)^2*(1.0+rt(i)^2/R(i)^2/((1.0+coeff_q*rt(i)^2/a^2)^2)));
Bphi(i)=B0/(1+rt(i)*cos(theta(i))/R0);
Btheta(i)=rt(i)*Bphi(i)/((1+8*rt(i)^2/a^2)*R(i));
end

% Calcul de la vitesse perpendiculaire 
vperp0=1;
vperp= zeros(N,1);
for i=1:N
vperp(i)=sqrt((vperp0^2*B(i))/B0);
end

% Résolution du système : x(1):rt, x(2):theta, x(3):phi, x(4):v//

dx= zeros(4,N);
x= zeros(4,N);
for i=1:N
dx(1,i) = -m(i)*R0*B0/R(i)^2/q(i)/B(i)^2*(x(4,i)^2+0.5*vperp(i)^2)*sin(x(2,i))+(Ett(i)*Bphi(i)-Ept(i)*Btheta(i))/B(i)^2;
dx(2,i) = x(4,i)/(1.0+coeff_q*x(1,i)^2/a^2)/R0-m(i)*B0*R0/R(i)^2/x(1,i)/q(i)/B(i)^2*(x(4,i)^2+0.5*vperp(i)^2)*cos(x(2,i))-(Ert(i)*Bphi(i)/rt(i));
dx(3,i)=x(4,i)/(R0+x(1,i)*cos(x(2,i)))+(Ert(i)*Btheta(i)/(R(i)*B(i)^2));    
dx(4,i)=-m(i).*vperp(i)^2/2.0/B(i)./m(i)./(1.0+coeff_q*x(1,i)^2/a^2)/(R(i)^3)*(x(1,i)*R0*B0*sin(x(2,i)))+(q(i)*sum(Ettot(i).*B(i))/(m(i)*B(i)));
end
end


Merci d'avance pour votre aide!