Equazioni differenziali a coefficienti non costanti

di il
473 risposte

473 Risposte - Pagina 10

  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    è ys che ha solo l'ultima colonna non nulla, non è un problema dell'assegnazione
    no
    questa è l'ultima ys:
    0 0 0 0,714565943272305 -0,401895682320436 2353163,47028039
    3,03661918715106e-07 -1,70789575357964e-07 1,00000026752371 0,714565943272305 -0,401895682320436 2353163,47028039


    l'ultimo elemento di T
    ha tutti zeri tranne l'ultima colonna che è uguale all'ultima colonna di ys
    0 0 0 0 0 2353163,47028039
    0 0 0 0 0 2353163,47028039
  • Re: Equazioni differenziali a coefficienti non costanti

    Ahhhh

    ho capito

    come puoi notare prima delle varie val(:,:,1) scrive così
    1.0e+006 *

    il che vuol dire che quelli non sono zeri...
    che se li moltipichi per 10^6 dentro hanno dei numeri...
    è solo una questione di visualizzazione
    Perché l'ultima colonna ha una scala diversa...
    ma i numeri ci sono
  • Re: Equazioni differenziali a coefficienti non costanti

    Prova, ad esempio, a guardarne uno per volta
    T(2,1,1)
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    Prova, ad esempio, a guardarne uno per volta
    T(2,1,1)

    perchè hai sempre ragione tu....
  • Re: Equazioni differenziali a coefficienti non costanti

    Lo so, lo so... è difficile da accettare...
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    Lo so, lo so... è difficile da accettare...

    comunque non le ammazza lo stesso
  • Re: Equazioni differenziali a coefficienti non costanti

    Il problema è sempre che non entra mai nell'if
    se fai min(T(end,1,:))
    è minore di r_pin1
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    il problema è sempre che non entra mai nell'if
    se fai min(T(end,1,:))
    è minore di r_pin1
    l'ho notato, quindi fa passare tutto...

    anche se gli chiedo il massimo, sono tutti numero piccolissimi...

    allora c'è qualcosa che non va prima
  • Re: Equazioni differenziali a coefficienti non costanti

    Ad occhio direi di sì...
  • Re: Equazioni differenziali a coefficienti non costanti

    Ieri sera, mentre ero a lavoro, ho capito cosa faccio di sbagliato...
    Io già che sono cretino, poi mi fanno fare 100 cose contemporaneamente e non capisco più nulla peggio del solito e faccio minchiate.

    Sistemo e mi fai un testing? grazie...

    P.S.: ho visto che con le strutture posso fare cose più belle che con le matrici.
  • Re: Equazioni differenziali a coefficienti non costanti

    Sicuramente puoi fare molte più cose... puoi, ad esempio, creare matrici di matrici o inserire in ciascuna cella un tipo di dato diverso (stringhe, numeri, matrici) e puoi anche dare un nome a ciascun campo, il che può risultare comodo nel momento in cui devi andare a recuperare i valori che ti servono.
  • Re: Equazioni differenziali a coefficienti non costanti

    La storia era che mi servono le componenti della velocità lungo i tre assi, non una velocità lungo un asse e due angoli per le altre due componenti, che è una cosa senza senso....

    Allora, 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=100;
    
    Ein=3;
    m_ion=1.00794*uma;
    
    q_over_m=q/m_ion;
    
    Etot_ion=zeros(n,1);
    betasquare_ion=zeros(n,1);
    v_ion=zeros(n,1);
    
    %%%%%%%%%%%%
    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.
    
    %T=zeros(2,6,n); % matrice tridimensionale dove immagazino le traiettore
    
    %calcolo la componente z della velocità
    for i=1:n              
             K_ion=normrnd(Ein,Ein/10);
             Etot_ion(i)=K_ion+((m_ion*uma1)/uma);  % ho tutto in MeV
             betasquare_ion(i)=1-(((m_ion*uma1)/(uma*Etot_ion(i)))^2);
             v_ion(i)=sqrt(betasquare_ion(i)*c^2);
             
             
             %spread angolare 
             A = -45; B = 45;
             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(fi)*sind(teta);
             
              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];
              
              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);
              
              %struttura
         Particella(i).traiettoria=ys(:,1:3);
         Particella(i).velocita=ys(:,4:6);
      
    end
    
    per ammazzare le particelle gli faccio fare un altro ciclo:
    
        for i=1:n
                if ( abs(Particella(i).traiettoria(end,1))>= r_pin1 ||  abs(Particella(i).traiettoria(end,1))>= r_pin1)
                    Particella(i)=[];
                    n=n-1;
                    i=n-1;
                end
        end
    
    e un poco le ammazza perchè la struttura Particella (senti che termine da uno che ne capisce...) esce dal primo ciclo che è 1x100
    ma ad un certo punto ho un errore nel secondo ciclo:
    ??? Index exceeds matrix dimensions.
    
    Error in ==> Prova at 99
                if ( abs(Particella(i).traiettoria(end,1))>= r_pin1 ||
                abs(Particella(i).traiettoria(end,1))>= r_pin1)
    
    e la struttura si è ridotta a 1x50 ma lui si incazza (come sempre) perchè l'indice diventa più grosso della matrice... io però ho copiato il codice tuo....(e questo dimostra che non ci capisco nulla)

    Mi dici anche se con le strutture faccio giusto?
  • Re: Equazioni differenziali a coefficienti non costanti

    Allora...
    l'errore è nell'aver fatto un ciclo a parte (il primo errore).
    Il secondo errore è nell'aver copiato male quello che ho scritto:
    i=n-1 -> i=i-1
    Il terzo errore è nell'aver aggiunto una riga a piacere che io non avevo scritto:
    n=n-1
    A parte tutto, leggendo il codice mi sembra che si possano semplificare notevolmente le cose.
    Però devi rinunciare al piacere di dover uccidere materialmente la particella...
    Fai così, crei le caselle della struttura Particella direttamente dentro l'if. Quindi se la condizione viene soddisfatta, la struttura si popola, altrimenti si va oltre.
    In questo caso è anche più corretto perchè invece che andare a modificare l'indice del ciclo usi un indice di supporto (k) che si incrementa solo se è verificata la condizione (lo devi inizializzare a 1 all'inizio del codice)
    if ( abs(ys(:,1:3))>= r_pin1 || abs(ys(:,1:3))>= r_pin1)
    Particella(k).traiettoria=ys(:,1:3);
    Particella(k).velocita=ys(:,4:6);
    k=k+1
    end
    Non l'ho provato ma dovrebbe funzionare. Matlab è più contento e le particelle che non soddisfano la condizione anche...
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    Allora...
    l'errore è nell'aver fatto un ciclo a parte (il primo errore).
    Il secondo errore è nell'aver copiato male quello che ho scritto:
    i=n-1 -> i=i-1
    Il terzo errore è nell'aver aggiunto una riga a piacere che io non avevo scritto:
    n=n-1
    Questa roba è uscita dopo vari tentativi...

    giug ha scritto:


    A parte tutto, leggendo il codice mi sembra che si possano semplificare notevolmente le cose.
    Però devi rinunciare al piacere di dover uccidere materialmente la particella...
    Fai così, crei le caselle della struttura Particella direttamente dentro l'if. Quindi se la condizione viene soddisfatta, la struttura si popola, altrimenti si va oltre.
    In questo caso è anche più corretto perchè invece che andare a modificare l'indice del ciclo usi un indice di supporto (k) che si incrementa solo se è verificata la condizione (lo devi inizializzare a 1 all'inizio del codice)
    if ( abs(ys(:,1:3))>= r_pin1 || abs(ys(:,1:3))>= r_pin1)
    Particella(k).traiettoria=ys(:,1:3);
    Particella(k).velocita=ys(:,4:6);
    k=k+1
    end
    Non l'ho provato ma dovrebbe funzionare. Matlab è più contento e le particelle che non soddisfano la condizione anche...
    Stavolta ho copia/incollato prima dell'end del ciclo for e lui mi dece:

    ??? Operands to the || and && operators must be convertible to
    logical scalar values.

    Error in ==> Prova at 86
    if ( abs(ys(:,1:3))>= r_pin1 || abs(ys(:,1:3))>= r_pin1)
  • Re: Equazioni differenziali a coefficienti non costanti

    Ho anche detto k=1; prima del for e corretto k=k+1 in k=k+1; perchè l'ho visto che si stava innervosendo...
Devi accedere o registrarti per scrivere nel forum
473 risposte