Equazioni differenziali a coefficienti non costanti

di il
473 risposte

473 Risposte - Pagina 16

  • Re: Equazioni differenziali a coefficienti non costanti

    Se non sbaglio si può fare l'eseguibile con il matlab compiler... non l'ho mai fatto...
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    Se non sbaglio si può fare l'eseguibile con il matlab compiler... non l'ho mai fatto...
    Poi lo faccio
  • Re: Equazioni differenziali a coefficienti non costanti

    Ho un problema con il plot finale...
    Voglio fare la storia dei doppi assi e gli faccio fare il plot così:
    
     %Contorno MCP
     Rmcp=0.02; %Raggio MCP
    Xmcporigine=linspace(0,Rmcp,1000); %crea un vettore riga di estremi 0-Rmpc costituito da 1000 elementi equidistanti
    Ymcporigine=sqrt(Rmcp^2-Xmcporigine.^2); %Equazione circonferenza MCP
    
    
    figure()
    hold all
    for i=1:p
        B_defl(i) = Part_trasp(i).traiettoria(end,2);
        E_defl(i) = Part_trasp(i).traiettoria(end,1);
        
        hplot = plot(B_defl(i), E_defl(i),'*','Markersize', 10);
    end
    
    hold on
    plot(Ymcporigine,Xmcporigine,'k-','MarkerEdgeColor','k','MarkerSize',1);  %Plot contorno MCP
    
    % xlabel('Magnetic deflection')
    % ylabel('Electric deflection')
     grid on
     
         %Label assi principali - It is necessary to give the label instructions after plot in order to avoid overlap
        xlabel(gca, 'Deflessione magnetica [m]'); % label lower x axis
        ylabel(gca,'Deflessione elettrica [m]');  %label left y axis
        
        %particles outside MCP radius won't appear in figure
        xlim([0, Rmcp])
        ylim([0, Rmcp])
        
    %     %%%% legenda assi principali
    %     l=legend(legendtext, 'Location', 'BestOutside'); 
    %     a=get(l, 'children');
    %     set(a(1:3:end),'MarkerSize', 10);
    % %     
    %   %%%% doppi Assi
        
        %set secondary x limit as the momentum of a H+ at distance equal to the MCP radius
        % Secondo Harres y=(q*B*LB*L)/sqrt(2mEkin) ==> mv=q*B*LB*L/y
        L=Deriva2;
        mv_max = (q*B*LB*L)/Rmcp;
        %mv_max = 1;
        
        %set secondary y limit as the energy of a H+ at distance equal to the MCP radius
        % Secondo Harres x=(q*E*Le*L)/(2Ekin) ==> Ekin=q*E*Le*L/2x
        Le = 0.07; %Estensione del C.E.  : 70 mm
        E = 100000; %campo TP.m
        Ekin_max = ((q*E*Le*L)/(2*Rmcp))/eV;
        %mv_max = 1;
        
        %Layout instruction
        set(gca,'Box','off');   % Turn off the box surrounding the whole axes
        %axesUnits=get(gca,'Units');
        axesPosition = get(gca,'Position');          %# Get the current axes position
        hNewAxes = axes('Position', axesPosition,...  %# Place a new axes on top...
            'Units', 'normalized',...
            'ActivePositionProperty', 'Position',... %prova anche OuterPosition
            'Color','none',...           %#   ... with no background color
            'YAxisLocation', 'right',...
            'XAxisLocation','top',...    %#   ... located on the top
            'Ylim', [0, Ekin_max],...
            'Xlim', [0, mv_max],...      %#   ... should define x axis scale (need to set xmax = mv_max)
            'Box','off');                %#   ... and no surrounding box
        
        xlabel(hNewAxes,'Momentum (H^+)');  %# Add a label to the top axis
        ylabel(hNewAxes,'Energy (H^+)');  %# Add a label to the right axis
        set(gca, 'XTickLabel', num2str(get(gca,'XTick')','%g'))
        set(gca, 'YTickLabel', num2str(get(gca,'YTick')','%g'))
    
    Praticamente il codice che avevamo fatto...
    però mi dà errore qua: hNewAxes = axes('Position', axesPosition,...

    dice che:

    MatLab ha scritto:


    ??? Error using ==> horzcat
    CAT arguments dimensions are not consistent.
  • Re: Equazioni differenziali a coefficienti non costanti

    Ekin_max e mv_max sono due numeri?
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    Ekin_max e mv_max sono due numeri?

    NO!


    Non la sopporto più sta cosa che hai sempre ragione tu...
    Io ci sbatto 3 ore, tu lo guardi un secondo e risolvi... sono solo invidioso...

    Comunque ho corretto, basta dargli la componente del campo, non tutto il campo :p
  • Re: Equazioni differenziali a coefficienti non costanti

    Ora fa anche il resize giusto... forse perchè ancora non c'è la legenda

    P.S.: ovviamente scherzavo, non è che te la prendi?
  • Re: Equazioni differenziali a coefficienti non costanti

    No no, tranquillo, ci sono abituata (ad avere ragione...)
    Comunque è questione di abitudine, in quel caso l'errore poteva essere solo lì... è l'unica parte della funzione che poteva avere errori di concatenazione perché è l'unica dove sono definiti dei vettori.
  • Re: Equazioni differenziali a coefficienti non costanti

    Ho fatto un run con centomila particelle e sostanzialmente con cambia molto il risultato finale.

    La domanda è questa: se nella parte dove fa il trasporto dentro lo spettrometro faccio un unico cilco anzicchè uno per ogni zona, che succede? dovrebbe risolvere le equazioni per una particella, trasportarla per tutto il percorso e poi ricominciare, giusto?
    Si incazza? a livello di programmazio è una soluzione migliore?
  • Re: Equazioni differenziali a coefficienti non costanti

    Se i cicli hanno gli stessi indici e per i cicli successivi non è necessario che quelli precedenti siano conclusi (ad esempio perché quelli precedenti devono costruire una matrice che serve "tutta" per quelli successivi), allora puoi mettere tutto insieme.
    Se in un ciclo ti crei la matrice delle traiettorie e nel successivo ti serve l'ultima traiettoria del ciclo precedente, allora non si può fare. Ma se ti servono solo i dati che vengono creati di volta in volta allora sì.
    Mi sono spiegata abbastanza male?
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    Se in un ciclo ti crei la matrice delle traiettorie e nel successivo ti serve l'ultima traiettoria del ciclo precedente, allora non si può fare. Ma se ti servono solo i dati che vengono creati di volta in volta allora sì.
    Mi sono spiegata abbastanza male?

    No, ti sei spiegata bene.
    Ogni volta che passo da una zona alla successiva mi servono gli ultimi dati della soluzione dell'equazione per darli come condizione iniziale per la successiva... quindi non lo posso fare.
  • Re: Equazioni differenziali a coefficienti non costanti

    Esatto
  • Re: Equazioni differenziali a coefficienti non costanti

    Eppure perchè dici così.
    Se io scrivo così
    
    for i=1:p   
        
        B=[0 0 0]'; 
        E=[0 0 0]';
          
        y0(i,:) = [Part_trasp(i).traiettoria(end,:),Part_trasp(i).velocita(end,:)];      
        %%%Risolvo il moto
        f = @(t,y1) [y1(4:6); (q_over_m).*cross(y1(4:6),B)+(q_over_m).*E]; 
                                                                                              
        options=odeset('RelTol',1e-7,'Events',@(t,y)Event_Stop_1(t,y,Deriva1));   
         [t,y1] = ode23t(f,tspan,y0(i,:),options);
         
        Part_trasp(i).traiettoria(end+1:end+size(y1,1),1:3)=y1(:,1:3);   
        Part_trasp(i).velocita(end+1:end+size(y1,1),1:3)=y1(:,4:6);  
    end
    %% 
    %Zona 2
    %calcolo la traiettoria nella parte in cui c'è solo il campo magnetico 
    for i=1:p
         
        B2=[0.003 0 0]';
        E2=[0 0 0]';
        
        y0(i,:) = [Part_trasp(i).traiettoria(end,:),Part_trasp(i).velocita(end,:)];  
         %%%Risolvo il moto
        f = @(t,y2) [y2(4:6); (q_over_m).*cross(y2(4:6),B2)+(q_over_m).*E2]; 
                                                  
        options=odeset('RelTol',1e-7,'Events',@(t,y)Event_Stop_2(t,y,inizioE));             
             
        [t,y2] = ode23t(f,tspan,y0(i,:),options);
        
        Part_trasp(i).traiettoria(end+1:end+size(y2,1),1:3)=y2(:,1:3);  
        Part_trasp(i).velocita(end+1:end+size(y2,1),1:3)=y2(:,4:6); 
    end
    
    Lui entra nel secondo ciclo, si crea la matrice delle condizioni iniziali (una riga per particella) e risolve.
    uso un solo ciclo lui che fa?
    risolve la prima equazione per i-esima particella, aggiorna la struttura Particella(i) e poi entra nella seconda zona e fa:
    y0(i,:) = y0(i,:) = [Part_trasp(i).traiettoria(end,:),Part_trasp(i).velocita(end,:)];
    ma solo per l'iesima particella, non per tutte, quindi dovrebbe creare un vettore riga per le condizioni iniziali che usa per risolvere l'altra equazione... o no?
  • Re: Equazioni differenziali a coefficienti non costanti

    Sì, è giusto, perché la i la usi solo per identificare la particella, non la usi per identificare gli elementi della traiettoria. E dato che l'i-esima particella non influenza la i+1, puoi mettere tutto insieme. ti basta cancellare qualche "for i=1:p" e qualche "end", il resto dovrebbe rimanere invariato.
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    sì, è giusto, perché la i la usi solo per identificare la particella, non la usi per identificare gli elementi della traiettoria. E dato che l'i-esima particella non influenza la i+1, puoi mettere tutto insieme. ti basta cancellare qualche "for i=1:p" e qualche "end", il resto dovrebbe rimanere invariato.

    EVVIVA!
    (per la prima volta ho ragione io... facciamo come Giotto e Cimabue? )
  • Re: Equazioni differenziali a coefficienti non costanti

    Va bene, ora puoi aprire un forum di matlab tutto tuo...
Devi accedere o registrarti per scrivere nel forum
473 risposte