Implementazione RK impliciti su MATLAB

di
cecioho
il
1 risposte
C
cecioho
Pochi posts
Iscritto da
dic, 2020
Messaggi:
1

Implementazione RK impliciti su MATLAB

Ciao, ho dei problemi nell'implementazione su MATLAB di un runge-kutta implicito generico con il metodo di Newton, nel senso che a me i codici sembrano giusti, li ho ricontrollati più e più volte ma non funziona, qualcuno riesce a darmi una mano? A è la matrice di Butcher, b e c i vettori del metodo

function y = rkimplis(f, Jf, y0, t0, t1, h, c, b, A, toll, nMax)

n = ceil((t1-t0)/h);

t = t0;
s = length(c);
d = length(y0);

k = zeros(s,d);
kk = zeros(1,s*d);

y = zeros(n+1,d);
y(1,:) = y0;

for l=1:n
    iter = 0;
    err = 100;
    
    while(err>toll && iter<nMax)
        iter = iter + 1;
        
        J = eye(s*d,s*d);
        F = zeros(1,s*d);
        
        for i=1:s
            J_block = Jf(t+c(i)*h,y(l,:)+h*A(i,:)*k);
            %for j=1:s
                %J((i-1)*d+1:i*d,(j-1)*d+1:j*d) = J((i-1)*d+1:i*d,(j-1)*d+1:j*d) - h*A(i,j)*J_block;
            %end
            
            J((i-1)*d+1:i*d,:) = J((i-1)*d+1:i*d,:) - h*kron(A(i,:),J_block);
            
        end
        
        for i=1:s
            F(1,(i-1)*d+1:i*d) = - kk((i-1)*d+1:i*d) + f(t+c(i)*h,y(l,:)+h*A(i,:)*k);
        end
        
        kk_temp = J\F';
        
        err = norm(kk_temp,'inf');
        
        kk = kk_temp' + kk;
        
        for i=1:s
            k(i,:) = kk((i-1)*d+1:i*d);
        end
    end
    
    y(l+1,:) = y(l,:) + h*b*k;
    
    if(t+h>t1)
        h = t1 - t;
        t = t1;
    else
        t = t + h;
    end
end



return

1 Risposte

  • F
    feddy
    Pochi posts
    Iscritto da
    lug, 2020
    Messaggi:
    29

    Re: Implementazione RK impliciti su MATLAB

    Ciao,

    premetto che non ho testato nulla, ma sto dando un'occhiata qua e là. Dal momento che non stai implementando un RK a passo variabile, sarebbe meglio fare un for che scorre ogni istante temporale, evitando

    cecioho ha scritto:

    if(t+h>t1)
        h = t1 - t;
        t = t1;
    else
        t = t + h;
    end
    Vedo che quando aggiorni la soluzione hai 
    h*b*k
    ma b e k sono vettori, sei sicuro che faccia quello che desideri? Cioè il prodotto scalare? Se le dimensioni non sono corrette, potrebbe ritornare una matrice.


    Per essere certi dell'implementazione con un ODEsolver di solito si inizia a testare con casi banali come y'=-y, y(0)=1. E poi se ne verifica l'ordine di convergenza. Cosa ottieni?
Devi accedere o registrarti per scrivere nel forum
1 risposte