Equazioni differenziali a coefficienti non costanti

di il
473 risposte

473 Risposte - Pagina 27

  • Re: Equazioni differenziali a coefficienti non costanti

    Diciamo che ho finito di scrivere lo script, ora lo trasformo in funzione.
    Puoi dare un'occhiata se c'è qualcosa che non va?
    
    A=-1; B=1;
    v=A+(B-A)*rand(200000,3); %cubo di lato 2 di numeri random uniformemente distribuiti di lato 2
    
    % figure()
    % plot3(v(:,1),v(:,2),v(:,3),'.b','MarkerSize',1)
    
    d = v(:,1).^2 + v(:,2).^2 + v(:,3).^2; %equzione della sfera
    posizioni=find(d<=1 & v(:,3)>=0 & v(:,1)~=0 & v(:,2)~=0 ); %cerco quei punti che stanno nella semisfera di raggio 1 e con z>=0 e x,y ~=0
     
    condizione=v(posizioni,:); %salvo i punti che stanno nella sfera
     
    % figure()
    %   plot3(v(:,1),v(:,2),v(:,3),'.y','MarkerSize',1)
    %   hold on
    % plot3(condizione(:,1),condizione(:,2),condizione(:,3),'.k','MarkerSize', 1)
    
    modulo_quadro=condizione(:,1).^2+condizione(:,2).^2+condizione(:,3).^2; %normalizzo i punti della semisfera in modo da ottenere una superfice
    modulo=sqrt(modulo_quadro);
    
    versore=condizione./repmat(modulo,1,3) ; % Calcolo i versori
                                             % repmat: Replicate and tile array. B = repmat(A,m,n) creates a large matrix B 
                                             % consisting of an m-by-n tiling of copies of A. The size of B is [size(A,1)*m, (size(A,2)*n].
                                             % Devo dividere le tre colonne di condizione per il vettore modulo, 
                                             % quindi faccio una matrice che contenga tre volte il vettore modulo per poter fare la divisione 
                                             % elemento per elemento. repmat(modulo,1,3) è uguale a scrivere [modulo modulo modulo]. 
    
    figure()
    
    plot3(versore(:,3),versore(:,2),versore(:,1), '.', 'MarkerSize',1)
    axis equal
    grid on
    xlabel('beam direction')
    ylabel('y axis')
    zlabel('x axis')
    
    z=0.95;                             % angular spread %%%INPUT
    n=100;                              % numero diparticelle %%%INPUT
    
    vers=find(versore(:,3)>=z);         % Selezioni i punti che stanno in una calotta sferica
    vers_select=versore(vers,:);        % creo il vettore con i versori da usare per le proiezioni della velocità
    
    % figure()
    % 
    % plot3(vers_select(:,3),vers_select(:,2),vers_select(:,1), '.', 'MarkerSize',1)
    % axis equal
    % grid on
    % xlabel('beam direction')
    % ylabel('y axis')
    % zlabel('x axis')
    
    r = randperm(size(vers_select,1)); %creates a vector of random number between 1 and size(vers_select)
    
    v=vers_select(r(1:n),:); % versori per il calcolo delle componenti della velocità
    
    figure()
    plot3(v(:,3),v(:,2),v(:,1), '.', 'MarkerSize',1)
    axis equal
    
    v_ion=2.8e4;      %Modulo delle velocità %%%INPUT
    v0=v.*v_ion;      %Matrice delle velocità iniziali %%%OUTPUT
    
    figure()
    plot3(v0(:,3),v0(:,2),v0(:,1), '.', 'MarkerSize',1)
    axis equal
    
    io dovrei dare in input n,z, e v_ion, quest'ultimo è un vettore, e ricevere in uscita v0 che è la matrice delle proiezioni di v_ion
    quindi l'intestazione, se non sbaglio è:

    function [v0] = nome_funzione(n,z, v_ion)
  • Re: Equazioni differenziali a coefficienti non costanti

    Ho un problema:

    se faccio un il vettore v_ion 100X1 e lo moltiplico per la matrice vers 100X3 lui si incazza.
    potrei fargli fare v_ion 100x3 con repmat(v_ion, 100,3) e andrebbe bene per adesso che v_ion è un numero.
    nel main che avevamo fatto però v_ion è un vettore, ad esempio 100x1, con righe tutte diverse quindi non posso usare repmat.
    cosa mi conviene fare? replicare la prima colonna tre volte? come gli dico di fare il prodotto di una riga di v_ion per le tre colonne della stessa riga di vers senza usare cicli?
  • Re: Equazioni differenziali a coefficienti non costanti

    Repmat serve per replicare, in generale delle matrici (quindi anche vettori) il numero di volte che vuoi nella direzione che vuoi.
    Se devi replicare un vettore 100x1 3 volte in modo che diventi un vettore 100x3 (quindi con 3 colonne uguali, basta che scrivi:
    repmat(v_ion,1,3)
    quindi gli stai dicendo di ripeterlo per una volta sola nel senso delle righe e 3 volte nel senso delle colonne.
  • Re: Equazioni differenziali a coefficienti non costanti

    Ho fatto un main di prova:
    
    z=0.95;
    n=100;
    v_ion_modulo=2.867e4;
    v_ion=repmat(v_ion_modulo,n,3);
    
    v0=feval(@God_Speed,n,z,v_ion);
    
    che chiama la funzione che ho sistemato così:
    
    function [v0] = God_Speed(n,z, v_ion)
    
    A=-1; B=1;
    v=A+(B-A)*rand(200000,3); %cubo di lato 2 di numeri random uniformemente distribuiti di lato 2
    
    % figure()
    % plot3(v(:,1),v(:,2),v(:,3),'.b','MarkerSize',1)
    
    d = v(:,1).^2 + v(:,2).^2 + v(:,3).^2; %equzione della sfera
    posizioni=find(d<=1 & v(:,3)>=0 & v(:,1)~=0 & v(:,2)~=0 ); %cerco quei punti che stanno nella semisfera di raggio 1 e con z>=0 e x,y ~=0
     
    condizione=v(posizioni,:); %salvo i punti che stanno nella sfera
     
    % figure()
    %   plot3(v(:,1),v(:,2),v(:,3),'.y','MarkerSize',1)
    %   hold on
    % plot3(condizione(:,1),condizione(:,2),condizione(:,3),'.k','MarkerSize', 1)
    
    modulo_quadro=condizione(:,1).^2+condizione(:,2).^2+condizione(:,3).^2; %normalizzo i punti della semisfera in modo da ottenere una superfice
    modulo=sqrt(modulo_quadro);
    
    versore=condizione./repmat(modulo,1,3) ; % Calcolo i versori
                                             % repmat: Replicate and tile array. B = repmat(A,m,n) creates a large matrix B 
                                             % consisting of an m-by-n tiling of copies of A. The size of B is [size(A,1)*m, (size(A,2)*n].
                                             % Devo dividere le tre colonne di condizione per il vettore modulo, 
                                             % quindi faccio una matrice che contenga tre volte il vettore modulo per poter fare la divisione 
                                             % elemento per elemento. repmat(modulo,1,3) è uguale a scrivere [modulo modulo modulo]. 
    
    figure()
    
    plot3(versore(:,3),versore(:,2),versore(:,1), '.', 'MarkerSize',1)
    axis equal
    grid on
    xlabel('beam direction')
    ylabel('y axis')
    zlabel('x axis')
    
    % z=0.95;                             % angular spread %%%INPUT
    % n=100;                              % numero diparticelle %%%INPUT
    
    vers=find(versore(:,3)>=z);         % Selezioni i punti che stanno in una calotta sferica
    vers_select=versore(vers,:);        % creo il vettore con i versori da usare per le proiezioni della velocità
    
    % figure()
    % 
    % plot3(vers_select(:,3),vers_select(:,2),vers_select(:,1), '.', 'MarkerSize',1)
    % axis equal
    % grid on
    % xlabel('beam direction')
    % ylabel('y axis')
    % zlabel('x axis')
    
    r = randperm(size(vers_select,1)); %creates a vector of random number between 1 and size(vers_select)
    
    vers=vers_select(r(1:n),:); % versori per il calcolo delle componenti della velocità
    
    figure()
    plot3(vers(:,3),vers(:,2),vers(:,1), '.', 'MarkerSize',1)
    axis equal
    
    %  v_ion_modulo=2.8e4;      %Modulo delle velocità %%%INPUT
    % v_ion=repmat(v_ion_modulo,n,3);
    v0=vers.*v_ion;      %Matrice delle velocità iniziali %%%OUTPUT
    
    figure()
    plot3(v0(:,3),v0(:,2),v0(:,1), '.', 'MarkerSize',1)
    axis equal
    
    e funziona, a meno che tu non vedi errori...

    Mi rimane da capire come sistemare la domanda di prima.
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    repmat serve per replicare, in generale delle matrici (quindi anche vettori) il numero di volte che vuoi nella direzione che vuoi.
    Se devi replicare un vettore 100x1 3 volte in modo che diventi un vettore 100x3 (quindi con 3 colonne uguali, basta che scrivi:
    repmat(v_ion,1,3)
    quindi gli stai dicendo di ripeterlo per una volta sola nel senso delle righe e 3 volte nel senso delle colonne.
    Non c'è verso di convincerlo a fare il prodotto matrice(100x3)*vettore(100x1)?
    ovviamente come dici tu funziona: faccio fare un vettore 100x1 al main di prova (come quello che poi dovrò passargli) e lo "matricizzo" nella funzione
  • Re: Equazioni differenziali a coefficienti non costanti

    No, non c'è verso... (o non lo so fare io)
    comunque se fai così, passandogli il vettore e usando repmat nella funzione rimane tutto abbastanza pulito e quasi non te ne accorgi.
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    No, non c'è verso... (o non lo so fare io)
    l'affermazione tra parentesi implica la prima
    comunque ho messo tutto insieme e non funziona nulla
    e la cosa peggiore è che lui nemmeno si incazza
  • Re: Equazioni differenziali a coefficienti non costanti

    Dopo una "lunga" ricerca posso affermare che il repmat è il metodo da utilizzare in questi casi.
    Se matlab non dà segni di malcontento... ti tocca lanciare riga per riga e vedere cosa succede.
  • Re: Equazioni differenziali a coefficienti non costanti

    Lo sto facendo... e poi è lento da morire...
  • Re: Equazioni differenziali a coefficienti non costanti

    Avrai fatto troppi cicli for non necessari...
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    Avrai fatto troppi cicli for non necessari...
    no, ti giuro, solo quelli del main che c'erano da prima
  • Re: Equazioni differenziali a coefficienti non costanti

    Allora bisogna capirlo, povero matlab... deve fare un sacco di calcoli...
  • Re: Equazioni differenziali a coefficienti non costanti

    giug ha scritto:


    Allora bisogna capirlo, povero matlab... deve fare un sacco di calcoli...
    fa solo questo e non capisce niente... povero io
  • Re: Equazioni differenziali a coefficienti non costanti

    Nel dubbio che i cicli del vecchio codice potessero essere un po' pesanti, ho dato un'occhiata...
    Puoi alleggerire di molto soprattutto i cicli interni.
    Ci sono i due cicli
    for s=1:stati
    ...
    for i=1:n
    ...
    end 
    end
    Quasi tutto il codice che sta all'interno del ciclo più interno lo puoi portare fuori.
    Ad esempio, la riga:
    betasquare_ion(i)=1-(((m_ions*uma1)/(uma*Etot_ion(i)))^2);
            
    non è necessario crearla elemento per elemento, al variare di i. Se la scrivi così:
    betasquare_ion=1-(((m_ions*uma1)/(uma.*Etot_ion)).^2);
    Ti crei la matrice al di fuori del ciclo e se ti serve all'interno indicizzi la riga. E la maggior parte neanche dovrai richiamarle all'interno del ciclo perché servono solamente per il calcolo successivo.
    Ricordati di mettere i puntini davanti alle moltiplicazioni/divisioni/quadrati in cui moltiplichi per un vettore.
    Così alleggerisci enormemente il codice, evitando di fargli fare un mare di calcoli.
    Lo stesso per quando crei gli angoli random. Ti crei direttamente una matrice rand della dimensione che ti serve e poi la devi solo indicizzare.
    Credo che l'unica (o quasi) riga di codice che deve rimanere nel ciclo interno è quella dell'ode23.
  • Re: Equazioni differenziali a coefficienti non costanti

    Altra cosa, sempre nel ciclo con s=1:stati... ma a cosa serve? mi sembra che nessuna variabile viene indicizzata con s. Questo vuol dire che lui ripete i calcoli tante volte quanti sono gli stati e poi memorizza solo i valori relativi all'ultimo (probabilmente c'è qualcosa che non vedo, o sto leggendo un codice vecchio...).
Devi accedere o registrarti per scrivere nel forum
473 risposte