Equazioni differenziali a coefficienti non costanti

di Anonimizzato9724 il
473 risposte

473 Risposte - Pagina 9

  • Non riesco a far funzionare il tuo codice perché mi manca l'eventStop
  • 1keenan ha scritto:


    l'errore è prima, quando assegno T(:,:,i)=ys

    nel workspace mi rimane l'ultimo ys, ma in T non c'è niente uguale a ys
    T è tutta fatta di zeri, quasi, per forza che non le ammazza, passano tutte...
  • Bene... risolto l'arcano...
  • Dovrei averti mandato script e funzione, ma ho fatto uncasino perchè nel frattempo sto facendoaltre 100 cose, dimmi se ti è arrivato tutto
  • Mi sembra di si...
    Comunque, tanto per fare compagnia a matlab ed essere un pò fastidiosa anche io... non mi piace come hai fatto l'if... devi mettere un || (OR) al posto dell'elseif dato che entrambe le condizioni portano allo stesso risultato T(:,:,i)=[];
    poi mi sembra che manchi l'aggiornamento dell'indice che se non lo metti ti salta le righe. Ma magari questa parte non mi è arrivata...
  • Non l'ho messo sicuro l'aggiornamento dell'indice, perchè non ho idea di cosa intendi...

    (che vergogna....)
  • 1keenan ha scritto:


    Non l'ho messo sicuro l'aggiornamento dell'indice, perchè non ho idea di cosa intendi...

    (che vergogna....)

    Non mi aiuti più ora che ho detto sta cosa, vero?
  • Riguardando con calma sembra che T la faccia tutta nulla tranne l'ultima colonna di ogni foglio che mette il valore giusto di ys...
    Allora scrivo T(:,:,i)=ys(:,:)
    e vedoche succede
  • Intendevo questo
    poi decrementi il contatore:
    i=i-1
    mi sembra difficile da non capire... quindi sicuramente l'hai fatto e non ho ricevuto la riga di codice, dispersa nella rete...
  • giug ha scritto:


    Intendevo questo
    poi decrementi il contatore:
    i=i-1
    mi sembra difficile da non capire... quindi sicuramente l'hai fatto e non ho ricevuto la riga di codice, dispersa nella rete...

    allora avevo capito giusto... non l'ho messo perchè hai detto che si incazza, quindi pensavo di fargli calcolare il contatore per il ciclo successivo in base alla "profondità" di T dopo che ha ammazzato le cose che dovrebbe ammazzare ma non ammazza... ma te lo avevo detto...
  • In ogni caso penso che è l'assegnamento di T che non va bene... assegna correttamente l'ultima colonna, le altre le fa tutte 0...
  • Ecco, non vedo dove gli fai calcolare il contatore....
    Comunque è proprio il fatto di andare a modificare il contatore che non gli piace, quindi tanto vale aggiungere i=i-1 e non complicarsi ulteriormente le cose.
    Questo il mio codice posticcio e mi sembra che funzioni bene:
    qui crea la matrice
    for i =1:5
    a(:,:,i)=rand(4);
    end
    qui cancella una riga:
    for i =1:5
        if i == 2
            a(:,:,i)=[];
            i=i-1;
        end
    end
    (L'eventstop non c'era nel messaggio).
  • 
    function [value,isterminal,direction]=Event_Stop_Sorgente(t,ys, Pinhole1)
    z=ys(3);                       %Assegno la posizione del fascio a z
    cond_z=z-Pinhole1;              %%% calcola la condizione di passaggio dallo zero
     
    value=cond_z;                  % funzione evento 
    
     
    isterminal=1;         %%% Stop the integration when the particle reach the magnets.
    direction=0;          %%% Stop the integration when the zero crossing occurs from every direction.
    
    questo è l'eventstop
  • È ys che ha solo l'ultima colonna non nulla, non è un problema dell'assegnazione
  • Ok, lasciamo perdere la cancellazione delle righe...
    se io gli faccio risolve le eq. differenziali per ogni i e glidico di salvare i dati in T lui mi fa una matrice che non contiene le soluzioni...
    usando questo ciclo:
    
    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);
             v_z=v_ion;
             
             %calcolo la componente x              
             v_x= unifrnd(-0.785,0.785, n,1); %angolo compreso da -45 e +45
                                  % returns an array of random numbers generated from the continuous uniform distributions 
                                  % with lower and upper endpoints specified by A and B, respectively. 
                                  % If A and B are arrays, R(i,j) is generated from the distribution specified
                                  % by the corresponding elements of A and B. If either A or B is a scalar, it is expanded
                                  % to the size of the other input. R = unifrnd(A,B,m,n,...) or R = unifrnd(A,B,[m,n,...]) 
                                  %returns an m-by-n-by-... array. If A and B are scalars, all elements of R are generated from 
                                  %the same distribution. If either A or B is an array, they must be m-by-n-by-... 
              
              %calcolo la componente x                                  
              v_y= unifrnd(-0.785,0.785, n,1);
             
              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); (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,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);
              
              T(:,:,i)=ys(:,:); %Matrice tridimensionale delle traiettore di tutte le particelle
    
    lui fa una T come quella in figura

    mentre le soluzioni che calcola non hanno tutti elementi nulli a parte l'ultimo.
    Allegati:
    9724_48f8b31ca23721529291553a81a070d4.jpg
    9724_48f8b31ca23721529291553a81a070d4.jpg
Devi accedere o registrarti per scrivere nel forum
473 risposte