Runge-Kutta pour EDO couplés

Question

Je vais construire une fonction d'octave qui peut résoudre N couplé équation différentielle du type:

dx/dt = F(x,y,…,z,t) 
dy/dt = G(x,y,…,z,t) 
dz/dt = H(x,y,…,z,t) 

Avec l'une de ces trois méthodes (Euler, Heun et Runge-Kutta -4).

Le code suivant correspond à la fonction:

function sol = coupled_ode(E, dfuns, steps, a, b, ini, method) 
    range = b-a; 
    h=range/steps; 
    rows = (range/h)+1; 
    columns = size(dfuns)(2)+1; 
    sol= zeros(abs(rows),columns); 
    heun=zeros(1,columns-1); 
    for i=1:abs(rows) 
    if i==1 
     sol(i,1)=a; 
    else 
     sol(i,1)=sol(i-1,1)+h;  
    end 
    for j=2:columns 
     if i==1 
     sol(i,j)=ini(j-1); 
     else 
     if strcmp("euler",method) 
      sol(i,j)=sol(i-1,j)+h*dfuns{j-1}(E, sol(i-1,1:end));  
     elseif strcmp("heun",method) 
      heun(j-1)=sol(i-1,j)+h*dfuns{j-1}(E, sol(i-1,1:end));   
     elseif strcmp("rk4",method) 
      k1=h*dfuns{j-1}(E, [sol(i-1,1), sol(i-1,2:end)]); 
      k2=h*dfuns{j-1}(E, [sol(i-1,1)+(0.5*h), sol(i-1,2:end)+(0.5*h*k1)]); 
      k3=h*dfuns{j-1}(E, [sol(i-1,1)+(0.5*h), sol(i-1,2:end)+(0.5*h*k2)]); 
      k4=h*dfuns{j-1}(E, [sol(i-1,1)+h, sol(i-1,2:end)+(h*k3)]); 
      sol(i,j)=sol(i-1,j)+((1/6)*(k1+(2*k2)+(2*k3)+k4));  
     end 
     end 
    end 
    if strcmp("heun",method) 
     if i~=1 
     for k=2:columns 
      sol(i,k)=sol(i-1,k)+(h/2)*((dfuns{k-1}(E, sol(i-1,1:end)))+(dfuns{k-1}(E, [sol(i,1),heun]))); 
     end 
     end 
    end  
    end 
end 

Lorsque j'utilise la fonction pour une seule équation différentielle ordinaire, la méthode RK4 est le meilleur que prévu, mais quand je courais le code pour un système couple de l'équation différentielle, RK4 est le pire, j'ai vérifié et vérifié et je ne sais pas ce que je fais mal.

Le code suivant est un exemple de comment appeler la fonction

F{1} = @(e, y) 0.6*y(3); 
F{2} = @(e, y) -0.6*y(3)+0.001407*y(4)*y(3); 
F{3} = @(e, y) -0.001407*y(4)*y(3); 

steps = 24; 

sol1 = coupled_ode(0,F,steps,0,24,[0 5 995],"euler"); 
sol2 = coupled_ode(0,F,steps,0,24,[0 5 995],"heun"); 
sol3 = coupled_ode(0,F,steps,0,24,[0 5 995],"rk4"); 

plot(sol1(:,1),sol1(:,4),sol2(:,1),sol2(:,4),sol3(:,1),sol3(:,4)); 
legend("Euler", "Heun", "RK4"); 

Answer 1

Attention: il y a un peu trop de l » dans le RK4 formulæ h:

k2 = h*dfuns{ [...] +(0.5*h*k1)]); 
k3 = h*dfuns{ [...] +(0.5*h*k2]); 

doivent être

k2 = h*dfuns{ [...] +(0.5*k1)]); 
k3 = h*dfuns{ [...] +(0.5*k2]); 

(dernières « s enlevés h).

Cependant, cela ne fait aucune différence pour l'exemple que vous avez fourni, depuis h=1 là.

Mais autre que ce petit bug, je ne pense pas que vous faites en fait quelque chose de mal.

Si je conspire la solution générée par le plus avancé, d'adaptation 4ᵗʰ/5ᵗʰ pour RK mis en œuvre ode45:

F{1} = @(e,y) +0.6*y(3); 
F{2} = @(e,y) -0.6*y(3) + 0.001407*y(4)*y(3); 
F{3} = @(e,y)   -0.001407*y(4)*y(3); 

tend = 24; 
steps = 24; 
y0 = [0 5 995]; 
plotN = 2; 

sol1 = coupled_ode(0,F, steps, 0,tend, y0, 'euler'); 
sol2 = coupled_ode(0,F, steps, 0,tend, y0, 'heun'); 
sol3 = coupled_ode(0,F, steps, 0,tend, y0, 'rk4'); 

figure(1), clf, hold on 
plot(sol1(:,1), sol1(:,plotN+1),... 
    sol2(:,1), sol2(:,plotN+1),... 
    sol3(:,1), sol3(:,plotN+1)); 

% New solution, generated by ODE45 
opts = odeset('AbsTol', 1e-12, 'RelTol', 1e-12); 

fcn = @(t,y) [F{1}(0,[0; y]) 
       F{2}(0,[0; y]) 
       F{3}(0,[0; y])]; 
[t,solN] = ode45(fcn, [0 tend], y0, opts);  
plot(t, solN(:,plotN)) 

legend('Euler', 'Heun', 'RK4', 'ODE45'); 
xlabel('t');  

Ensuite, nous avons quelque chose de plus crédible à comparer.

Maintenant, RK4 plaine et simple, effectue en effet terriblement pour ce cas isolé:

Cependant, si je feuillette simplement les signes du dernier terme dans les deux dernières fonctions:

%      ± 
F{2} = @(e,y) +0.6*y(3) - 0.001407*y(4)*y(3); 
F{3} = @(e,y)   +0.001407*y(4)*y(3); 

Ensuite, nous obtenons ceci:

La principale raison pour laquelle RK4 fonctionne mal pour votre cas est en raison de la taille du pas. Le RK4/5 adaptatif (avec une tolérance fixée à 1 au lieu de 1e-12 comme ci-dessus) produit un δt moyen = 0,15. Cela signifie que l'analyse d'erreur de base a indiqué que pour ce problème particulier, h = 0.15 est la plus grande étape que vous pouvez prendre sans introduire d'erreur inacceptable.

Mais vous prenez h = 1, ce qui donne alors en effet une grande erreur accumulée.Le fait que Heun et Euler fonctionnent si bien pour votre cas est, bien, chance simple, comme démontré par l'exemple d'inversion de signe ci-dessus.

Bienvenue dans le monde des mathématiques numériques - il n'y a jamais une méthode qui convient mieux à tous les problèmes en toutes circonstances :)