Equazioni differenziali a coefficienti non costanti

di il
473 risposte

473 Risposte - Pagina 17

  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    Va bene, ora puoi aprire un forum di matlab tutto tuo...
    Ma nemmeno per sogno, significa imparare un minimo di php o quello che si usa per il web... io non la conosco sta roba e non la voglio conoscere.
    (se perdo la mia ignoranza poi chi me la ritorna?)

    Ho fatto un solo ciclo per la parte di trasporto, quella prima no, perchè perdo particelle e là credo che è necessario fare 2 cicli separati, il risultato grafico è bene o male lo stesso, certo quale badduzza è sparata un pochino più in là, ma credo dipenda dalla selezione random di angolo e energia.
  • Re: Equazioni differenziali a coefficienti non costanti

    Dovrebbe darti il risultato esattamente equivalente alla situazione con i cicli separati (a meno della casualità del processo).
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    Dovrebbe darti il risultato esattamente equivalente alla situazione con i cicli separati (a meno della casualità del processo).
    era quello che volevo dire, ma tu parli meglio :p
  • Re: Equazioni differenziali a coefficienti non costanti

    Era la badduzza che mi lasciava perplessa...
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    era la badduzza che mi lasciava perplessa...
    badduzza = pallina
    che sarebbe la particella
  • Re: Equazioni differenziali a coefficienti non costanti

    Sì sì, immaginavo... anche perché è da 245 post che non parliamo d'altro...
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    ... è da 245 post che non parliamo d'altro...
    Non è vero, al più 180
  • Re: Equazioni differenziali a coefficienti non costanti

    180 sono solo i tuoi...
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    180 sono solo i tuoi...
    mi ero dimenticato di considerare anche te nella discussione... :D

    vabbè sto facendo la storia degli stati di carica e la domanda è questa:

    Nella struttura particella posso inserire un campo .stato che a seconda del valore gli associa la corretta energia e massa?
    Penso di sì, ma suppongo che questo voglia dire che devo dare anche altri due campo, cioè energia e massa, giusto?
  • Re: Equazioni differenziali a coefficienti non costanti

    Direi che il discorso fila...
  • Re: Equazioni differenziali a coefficienti non costanti

    Quindi il primo ciclo dello script dovrebbe diventare
    
    for s=1:stati
      Ein_s=j*Ein;
           m_ionj=m_ion-(j*m_e);
    
    
       for i=1:n              
             K_ion=normrnd(Ein_s,Ein_s/10);
             Etot_ion(i)=K_ion+((m_ionj*uma1)/uma);  % ho tutto in MeV
             betasquare_ion(i)=1-(((m_ionj*uma1)/(uma*Etot_ion(i)))^2);
             v_ion(i)=sqrt(betasquare_ion(i)*c^2);
                      
             %angular spread  
             %A = -0.0575; B = 0.0575;
             
             A=-0.08; B=0.08;
             %A=0; B=0;
             angolo = A + (B-A) * rand(1); %random number in [A,B]
             
             teta=angolo; %angolo nel piano zy
             fi=angolo; %angolo con l'asse x
                      
             %%%%Componente z della velocità
             v_z= v_ion*cosd(teta)*cosd(fi);
              
             %calcolo la componente x              
             v_x= v_ion*sind(teta);
             
              %calcolo la componente y 
              v_y=v_ion*cosd(teta)*sind(fi);
             
              v0=[v_x, v_y, v_z]; %initial velocity
              r0=zeros(n,3);      %initial speed
              
              y_sorgente=[r0,v0];  %vector of initial condition
              
              tspan=[0,100];
                        
              O=[0 0 0]';
              %%%Risolvo il moto
              f = @(t,ys) [ys(4:6); O];
              % the expression [y(4:6); O]
              % combines two vectors of length 3 to make a vector of
              % length 6(as long as y is a column vector).
              
              options=odeset('RelTol',1e-7,'Events',@(t,ys)Event_Stop_Sorgente(t,ys,Pinhole1));   % opzione per risoluzione equazione differenziale:
              %set precisione calcolo---% Sintassi per utilizzare Event_Stop con Pinhole1 come variabile in input
              [t,ys] = ode23t(f,tspan,y_sorgente(i,:),options);
                        
              %decido quali particelle passano il pimo diaframma
              
         if (abs(ys(end,1)) <=r_pin1||abs(ys(end,2))<=r_pin1)
             
            % T(:,:,k)=ys;
                 
                 Particella(k).stato_di_carica=s;
                 Particella(k).massa=m_ionj;
                 Particella(k).q_over_m= s/m_ionj;
                 Particella(k).energia_iniziale=Ein_s;
                 Particella(k).traiettoria=ys(:,1:3);
                 Particella(k).velocita=ys(:,4:6);
                 k=k+1;
          end
      end      
    end
    
    lo provi e mi fai sapere se si incazza? io ho troppa paura
  • Re: Equazioni differenziali a coefficienti non costanti

    Ho sbagliato.
    deve essere:

    Particella(k).stato_di_carica=s*q;
    Particella(k).q_over_m= (s*q)/m_ionj;
  • Re: Equazioni differenziali a coefficienti non costanti

    Perchè escono fuori numeri complessi?

    L'ho capito forse: ho fatto un bordello con l'indice s che poi è diventato j... e questo l'ho corretto...
  • Re: Equazioni differenziali a coefficienti non costanti

    Dovrebbe essere così:
    
    clear all
    close all
    
    uma=1.66e-27; %[Kg], unità di massa atomica
    q=1.602e-19; %[C], carica elementare
    
    m_e=5.485799e-4*uma; %[kg], massa dell'elettrone
    eV=1.6021765e-19; %[J], fattore di conversione eV-->J
    c= 2.99792458e8;  %[m/s], velocità della luce
    %KB=(8.617332*10^(-5))*eV; %[eV/K], costante di Boltzmann
    uma1=931.494028; %[MeV/c^2], fattore di conversione uma-->MeV/c^2 (massa)
    KB=8.617e-5; %[eV*K^-1], costante di Boltzmann
    KBT=1076*KB;
    
    n=500; %numero diparticelle
    
    Ein=0.3;
    m_ion=1.00794*uma;
    
    %q_over_m=q/m_ion;
    
    stati= 3; %numero stati di carica
    
    Etot_ion=zeros(n,1);
    betasquare_ion=zeros(n,1);
    v_ion=zeros(n,1);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%
    %%
    %COLLIMAZIONE FASCIO
    %
    %%
    %%
    %Passaggio dal primo pinhole
    %%%%%%%%%%%%
    Pinhole1=1;
    DiametroPinhole1=0.001;
    r_pin1=DiametroPinhole1/2;
    
    %Have a for loop over a bunch of particles. For each particle, set its 
    %velocity (speed and angle) vector using the rand() function. Then use an 
    %"if" statement to decide whether the particle passed through or not and 
    %count those that do.
    
    k=1; %indice di supporto per selezionare le particelle nel ciclo interno if
    
    for s=1:stati
      Ein_s=s*Ein;
           m_ions=m_ion-(s*m_e);
    
       for i=1:n              
             K_ion=normrnd(Ein_s,Ein_s/10);
             Etot_ion(i)=K_ion+((m_ions*uma1)/uma);  % ho tutto in MeV
             betasquare_ion(i)=1-(((m_ions*uma1)/(uma*Etot_ion(i)))^2);
             v_ion(i)=sqrt(betasquare_ion(i)*c^2);
                      
             %angular spread  
             %A = -0.0575; B = 0.0575;
             
             A=-0.08; B=0.08;
             %A=0; B=0;
             angolo = A + (B-A) * rand(1); %random number in [A,B]
             
             teta=angolo; %angolo nel piano zy
             fi=angolo; %angolo con l'asse x
                      
             %%%%Componente z della velocità
             v_z= v_ion*cosd(teta)*cosd(fi);
              
             %calcolo la componente x              
             v_x= v_ion*sind(teta);
             
              %calcolo la componente y 
              v_y=v_ion*cosd(teta)*sind(fi);
             
              v0=[v_x, v_y, v_z]; %initial velocity
              r0=zeros(n,3);      %initial position
              
              y_sorgente=[r0,v0];  %vector of initial condition
              
              tspan=[0,100];
                        
              O=[0 0 0]';
              %%%Risolvo il moto
              f = @(t,ys) [ys(4:6); O];
              % the expression [y(4:6); O]
              % combines two vectors of length 3 to make a vector of
              % length 6(as long as y is a column vector).
              
              options=odeset('RelTol',1e-7,'Events',@(t,ys)Event_Stop_Sorgente(t,ys,Pinhole1));   % opzione per risoluzione equazione differenziale:
              %set precisione calcolo---% Sintassi per utilizzare Event_Stop con Pinhole1 come variabile in input
              [t,ys] = ode23t(f,tspan,y_sorgente(i,:),options);
                        
              %decido quali particelle passano il pimo diaframma
              
         if (abs(ys(end,1)) <=r_pin1||abs(ys(end,2))<=r_pin1)
             
            % T(:,:,k)=ys;
                 
                 Particella(k).stato_di_carica=s;
                 Particella(k).massa=m_ions;
                 Particella(k).q_over_m= s/m_ions;
                 Particella(k).energia_iniziale=Ein_s;
                 Particella(k).traiettoria=ys(:,1:3);
                 Particella(k).velocita=ys(:,4:6);
                 k=k+1;
          end
      end
               
    end
    
    solo che credo che n non è più il numero di particelle, ma il numero di particelle per stato di carica.
    quindi se n=500 e stati=3 le particelle in totale sono 1500, giusto?
  • Re: Equazioni differenziali a coefficienti non costanti

    Allora, il primo pezzo l'ho fatto e sembra funzionare.
    Ho fatto così:
    
    clear all
    close all
    
    uma=1.66e-27; %[Kg], unità di massa atomica
    q=1.602e-19; %[C], carica elementare
    
    m_e=5.485799e-4*uma; %[kg], massa dell'elettrone
    eV=1.6021765e-19; %[J], fattore di conversione eV-->J
    c= 2.99792458e8;  %[m/s], velocità della luce
    %KB=(8.617332*10^(-5))*eV; %[eV/K], costante di Boltzmann
    uma1=931.494028; %[MeV/c^2], fattore di conversione uma-->MeV/c^2 (massa)
    KB=8.617e-5; %[eV*K^-1], costante di Boltzmann
    KBT=1076*KB;
    
    n=500; %numero di particelle per stato di carica
    
    Ein=0.3;
    m_ion=1.00794*uma;
    
    %q_over_m=q/m_ion;
    
    stati= 3; %numero stati di carica
    
    Etot_ion=zeros(n,1);
    betasquare_ion=zeros(n,1);
    v_ion=zeros(n,1);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%
    %%
    %COLLIMAZIONE FASCIO
    %
    %%
    %%
    %Passaggio dal primo pinhole
    %%%%%%%%%%%%
    Pinhole1=1;
    DiametroPinhole1=0.001;
    r_pin1=DiametroPinhole1/2;
    
    %Have a for loop over a bunch of particles. For each particle, set its 
    %velocity (speed and angle) vector using the rand() function. Then use an 
    %"if" statement to decide whether the particle passed through or not and 
    %count those that do.
    
    k=1; %indice di supporto per selezionare le particelle nel ciclo interno if
    
    for s=1:stati
      Ein_s=s*Ein;
           m_ions=m_ion-(s*m_e);
    
       for i=1:n              
             K_ion=normrnd(Ein_s,Ein_s/10);
             Etot_ion(i)=K_ion+((m_ions*uma1)/uma);  % ho tutto in MeV
             betasquare_ion(i)=1-(((m_ions*uma1)/(uma*Etot_ion(i)))^2);
             v_ion(i)=sqrt(betasquare_ion(i)*c^2);
                      
             %angular spread  
             %A = -0.0575; B = 0.0575;
             
             A=-0.08; B=0.08;
             %A=0; B=0;
             angolo = A + (B-A) * rand(1); %random number in [A,B]
             
             teta=angolo; %angolo nel piano zy
             fi=angolo; %angolo con l'asse x
                      
             %%%%Componente z della velocità
             v_z= v_ion*cosd(teta)*cosd(fi);
              
             %calcolo la componente x              
             v_x= v_ion*sind(teta);
             
              %calcolo la componente y 
              v_y=v_ion*cosd(teta)*sind(fi);
             
              v0=[v_x, v_y, v_z]; %initial velocity
              r0=zeros(n,3);      %initial position
              
              y_sorgente=[r0,v0];  %vector of initial condition
              
              tspan=[0,100];
                        
              O=[0 0 0]';
              %%%Risolvo il moto
              f = @(t,ys) [ys(4:6); O];
              % the expression [y(4:6); O]
              % combines two vectors of length 3 to make a vector of
              % length 6(as long as y is a column vector).
              
              options=odeset('RelTol',1e-7,'Events',@(t,ys)Event_Stop_Sorgente(t,ys,Pinhole1));   % opzione per risoluzione equazione differenziale:
              %set precisione calcolo---% Sintassi per utilizzare Event_Stop con Pinhole1 come variabile in input
              [t,ys] = ode23t(f,tspan,y_sorgente(i,:),options);
                        
              %decido quali particelle passano il pimo diaframma
              
         if (abs(ys(end,1)) <=r_pin1||abs(ys(end,2))<=r_pin1)
             
            % T(:,:,k)=ys;
                 
                 Particella(k).stato_di_carica=s;
                 Particella(k).massa=m_ions;
                 Particella(k).q_over_m= s/m_ions;
                 Particella(k).energia_iniziale=Ein_s;
                 Particella(k).traiettoria=ys(:,1:3);
                 Particella(k).velocita=ys(:,4:6);
                 k=k+1;
          end
      end
               
    end
    %%
    %Passaggio dal secondo pinhole
    Pinhole2=0.1;
    
    Distanzap1p2=Pinhole1+Pinhole2;
    DiametroPinhole2=0.0005;
    r_pin2=DiametroPinhole2/2;
    
    m=k-1; %numero di particelle che passano
    
    y_pin2=zeros(m,6);  %%%matrice condizione iniziale per il passaggio dal secondo pinhole
    % for i=1:m
    %     
    % %y_pin2(i,:)=[Particella(i).traiettoria(end,:),Particella(i).velocita(end,:)];  %vector of initial condition
    % %y_pin2(i,:)=T(end,:,i);
    % end
    
    j=1; %
    for i=1:m
        y_pin2(i,:)=[Particella(i).traiettoria(end,:),Particella(i).velocita(end,:)];  %vector of initial condition
              tspan=[0,100];
                        
              %O=[0 0 0]';
              %%%Risolvo il moto
              f = @(t,yp) [yp(4:6); O];
              % the expression [y(4:6); O]
              % combines two vectors of length 3 to make a vector of
              % length 6(as long as y is a column vector).
              
              options2=odeset('RelTol',1e-7,'Events',@(t,yp)Event_Stop_Coll(t,yp,Distanzap1p2));   % opzione per risoluzione equazione differenziale:
              %set precisione calcolo---% Sintassi per utilizzare Event_Stop con Pinhole1 come variabile in input
              [t,yp] = ode23t(f,tspan,y_pin2(i,:),options2);
                        
              %decido quali particelle passano il pimo diaframma
              
         if (abs(yp(end,1)) <=r_pin2 || abs(yp(end,2))<=r_pin2)
             
            % T2(:,:,j)=yp;
                 Part_collimata(j).stato_di_carica =Particella(j).stato_di_carica;
                 Part_collimata(j).massa = Particella(j).massa;
                 Part_collimata(j).q_over_m= Particella(j).q_over_m;
                 Part_collimata(j).traiettoria=yp(:,1:3);
                 Part_collimata(j).velocita=yp(:,4:6);
                 j=j+1;
          end
    end
    
    %number of particles entering the spectrometer
    p=j-1;
    figure('color', 'white')
    hold all
    
    
    for j=1:m
        plot3(Particella(j).traiettoria(:,3),Particella(j).traiettoria(:,2), Particella(j).traiettoria(:,1))
        %F(j)= getframe;
    end
    for i=1:p
       
        
        plot3(Part_collimata(i).traiettoria(:,3),Part_collimata(i).traiettoria(:,2), Part_collimata(i).traiettoria(:,1))
        
        Part_trasp(i) = Part_collimata(i);  %particelle da trasportare
        
    end
    grid on
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%
    %%
    %%TRASPORTO DENTRO LO SPETTROMETRO
    %
    %%
    %%
    %Distanze in metri
    
    P2_B = 0.091; %distanza secondo pinhole campo magnetico (Deriva)
    Deriva1=Distanzap1p2+P2_B; %distanza totale percorsa dalla particella prima del settore di deflessione
    
    LB=0.15;   %Lunghezza campo magnetico
    
    P2_E = 0.191; %distanza secondo pinhole campo elettrico (solo campo magnetico)
    inizioE=Distanzap1p2+P2_E; %distanza totale percorsa dalla particella prima dell'ingresso nel campo elettrico
    
    P2_fineB=0.241; %distanza secondo pinhole fine campo magnetico (zona a due campi)
    fineB= Distanzap1p2+P2_fineB; %distanza totale percorsa dalla particella fino all'uscita dal campo magnetico
    
    LE=0.07; %Lunghezza campo elettrico
    
    P2_fineE=0.261; %distanza secondo pinhole fine campo elettrico (solo campo elettrico)
    fineE=Distanzap1p2+P2_fineE; %distanza totale percorsa dalla particella fino all'uscita dal campo elettrico
    
    Deriva2=0.393;  %Distanza spettrometro rivelatore
    Dist_Spettr_riv=Distanzap1p2+Deriva2; %distanza totale percorsa dalla particella fino al rivelatore
    
    %%
    %Zona 1
    %Calcolo traiettoria nel tratto di deriva tra collimatori e spettrometro
    y0=zeros(p,6); %inizializzo la matrice delle condizioni iniziali
    a=1;
    for i=1:p   
        %Campi nulli
        B=[0 0 0]'; 
        E=[0 0 0]';
          
        y0(i,:) = [Part_trasp(i).traiettoria(end,:),Part_trasp(i).velocita(end,:)];   %vettore delle condizini iniziali    
        
        %%%Risolvo il moto
        f = @(t,y1) [y1(4:6); (Part_trasp(i).q_over_m).*cross(y1(4:6),B)+(Part_trasp(i).q_over_m).*E]; 
                                                 % the expression [y(4:6); (q/m)*cross(y(4:6),B)]
                                                 % combines two vectors of length 3 to make a vector of 
                                                 % length 6(as long as y is a column vector).
                                                 
        options=odeset('RelTol',1e-7,'Events',@(t,y)Event_Stop_1(t,y,Deriva1));   % opzione per risoluzione equazione differenziale:
                    %set precisione calcolo---% Sintassi per utilizzare Event_Stop con Deriva1 come variabile in input
        [t,y1] = ode23t(f,tspan,y0(i,:),options);
           
        %traiettoria= y1; % aggiorno la traiettoria
        
        Part_trasp(i).traiettoria(end+1:end+size(y1,1),1:3)=y1(:,1:3);  %aggiungo alla traiettoria le nuove coordinate partendo dall'utlima riga
        Part_trasp(i).velocita(end+1:end+size(y1,1),1:3)=y1(:,4:6);  %aggiungo alla veloctà le nuove coordinate partendo dall'utlima riga
        
    end
    
    %PLOT
    figure('color', 'white')
    hold all
    
    for i=1:p
       
        plot3(Part_trasp(i).traiettoria(:,3),Part_trasp(i).traiettoria(:,2), Part_trasp(i).traiettoria(:,1))
        
        %Part_trasp(i) = Part_collimata(i);  %particelle da trasportare
        
    end
    grid on
    
    Poi nella parte successiva si blocca: dice "busy" e non da segni di vita. pensavo si bloccasse nel ciclo ma anche se dico:
    
    i=1;
     B2=[0.003 0 0]';
        E2=[0 0 0]';
        
        y0(i,:) = [Part_trasp(i).traiettoria(end,:),Part_trasp(i).velocita(end,:)];   %vettore delle condizini iniziali
        
         %%%Risolvo il moto
        f = @(t,y2) [y2(4:6); (Part_trasp(i).q_over_m).*cross(y2(4:6),B2)+(Part_trasp(i).q_over_m).*E2]; 
                                                 % the expression [y(4:6); (q/m)*cross(y(4:6),B)]
                                                 % combines two vectors of length 3 to make a vector of 
                                                 % length 6(as long as y is a column vector).
        
        options=odeset('RelTol',1e-7,'Events',@(t,y)Event_Stop_2(t,y,inizioE));            % opzione per risoluzione equazione differenziale, set precisione calcolo                                         
        [t,y2] = ode23t(f,tspan,y0(i,:),options);
        
        %traiettoria= [y1;y2];  % aggiorno la traiettoria
        Part_trasp(i).traiettoria(end+1:end+size(y2,1),1:3)=y2(:,1:3);  %aggiungo alla traiettoria le nuove coordinate partendo dall'utlima riga
        Part_trasp(i).velocita(end+1:end+size(y2,1),1:3)=y2(:,4:6);  %aggiungo alla veloctà le nuove coordinate partendo dall'utlima riga
    
    lui non la smette più di fare conti, e non capisco perchè
Devi accedere o registrarti per scrivere nel forum
473 risposte