Absoluter Fehler von ODE45 und Runge-Kutta-Methoden im Vergleich zur analytischen Lösung

Question

Ich würde mich freuen, wenn jemand mit folgendem Problem helfen kann. Ich habe folgende ODE:

dr/dt = 4*exp(0.8*t) - 0.5*r ,r(0)=2, t[0,1]  (1) 

I (1) auf zwei verschiedene Arten gelöst haben. Mittels der Runge-Kutta Methode (4. Ordnung) und mittels ode45 in Matlab. Ich verglich haben die beiden Ergebnisse mit der analytischen Lösung, die gegeben ist durch:

r(t) = 4/1.3 (exp(0.8*t) - exp(-0.5*t)) + 2*exp(-0.5*t) 

Als ich den absoluten Fehler jedes Verfahren in Bezug auf die exakte Lösung plotten, erhalte ich die folgende:

Für RK -Methode, mein Code:

h=1/50;            
x = 0:h:1;           
y = zeros(1,length(x)); 
y(1) = 2;  
F_xy = @(t,r) 4.*exp(0.8*t) - 0.5*r;     
for i=1:(length(x)-1)        
    k_1 = F_xy(x(i),y(i)); 
    k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1); 
    k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2)); 
    k_4 = F_xy((x(i)+h),(y(i)+k_3*h)); 
    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h; % main equation 
end 

Und für ode45:

tspan = 0:1/50:1; 
x0 = 2; 
f = @(t,r) 4.*exp(0.8*t) - 0.5*r; 
[tid, y_ode45] = ode45(f,tspan,x0); 

Meine Frage ist, warum habe ich Schwingungen, wenn ich ode45 verwenden? (Ich beziehe mich auf den absoluten Fehler). Beide Lösungen sind genau (1e-9), aber was passiert mit ode45 in diesem Fall?

Wenn ich den absoluten Fehler für die RK-Methode berechne, warum sieht es schöner aus?

Answer 1

Ihre RK4-Funktion nimmt feste Schritte vor, die viel kleiner sind als diejenigen, die ode45 dauert. Was Sie wirklich sehen, ist der Fehler aufgrund polynomial interpolation, der verwendet wird, um die Punkte zwischen den wahren Schritten zu produzieren, die ode45 dauert. Dies wird oft als "dichte Ausgabe" bezeichnet (siehe Hairer & Ostermann 1990).

Wenn Sie einen TSPAN Vektor mit mehr als zwei Elementen angeben, geben Sie Matlab's ODE suite solvers feste Schrittgrößenausgabe aus. Dies bedeutet nicht, dass sie tatsächlich eine feste Schrittgröße verwenden oder dass sie die in Ihrer TSPAN angegebenen Schrittgrößen verwenden. Sie können die tatsächlichen Schrittgrößen verwendet sehen und immer noch die gewünschte feste Schrittgröße Ausgabe erhalten, indem ode45 Ausgang eine Struktur aufweisen und mit deval:

sol = ode45(f,tspan,x0); 
diff(sol.x) % Actual step sizes used 
y_ode45 = deval(sol,tspan); 

Sie, dass nach einem ersten Schritt von 0.02 sehen werden, weil Ihre ODE ist einfach konvergiert es zu 0.1 für die nachfolgenden Schritte. Die Standardtoleranzen in Kombination mit der standardmäßigen maximalen Schrittgrößenbeschränkung (ein Zehntel des Integrationsintervalls) bestimmen dies.Schauen wir uns die wahren Schritte, um den Fehler Grundstück:

exactsol = @(t)(4/1.3)*(exp(0.8*t)-exp(-0.5*t))+2*exp(-0.5*t); 
abs_err_ode45 = abs(exactsol(tspan)-y_ode45); 
abs_err_ode45_true = abs(exactsol(sol.x)-sol.y); 
abs_err_rk4 = abs(exactsol(tspan)-y); 
figure; 
plot(tspan,abs_err_ode45,'b',sol.x,abs_err_ode45_true,'k.',tspan,abs_err_rk4,'r--') 
legend('ODE45','ODE45 (True Steps)','RK4',2) 

Wie Sie die Fehler an den wahren Schritte sehen können, langsamer als der Fehler für RK4 wächst (ode45 ist effektiv ein höherer Ordnung Verfahren als RK4 also würdest du das erwarten). Der Fehler wächst aufgrund der Interpolation zwischen den Integrationspunkten. Wenn Sie dies einschränken möchten, sollten Sie die Toleranzen oder andere Optionen über odeset anpassen.

Wenn Sie ode45 zu verwenden, um einen Schritt von 1/50 zwingen wollten, dies zu tun (funktioniert, weil Ihre ODE einfach ist):

opts = odeset('MaxStep',1/50,'InitialStep',1/50); 
sol = ode45(f,tspan,x0,opts); 
diff(sol.x) 
y_ode45 = deval(sol,tspan); 

Für ein weiteres Experiment, versuchen Sie das Integrationsintervall zu vergrößern, um t = 10 zu integrieren out könnte sein. Sie werden viel interessantes Verhalten in dem Fehler sehen (die Darstellung des relativen Fehlers ist hier nützlich). Kannst du das erklären? Können Sie ode45 und odeset verwenden, um Ergebnisse zu erzielen, die sich gut verhalten? Die Integration exponentieller Funktionen über große Intervalle mit adaptiven Schrittmethoden ist eine Herausforderung und ode45 ist nicht unbedingt das beste Werkzeug für den Job. Es gibt alternatives jedoch, aber sie erfordern möglicherweise einige Programmierung.

Answer 2

ode45 ist gekoppelt rk4-rk5. Persönlich denke ich, dass der ODE45-Fehler schöner ist. Beachten Sie, dass es begrenzt bleibt. Die ode4 wird korrigiert, wenn die Fehlergröße zu groß wird und der minimale Fehler pro Zyklus etwa 1e-10 ist. Das rk4 "rennt weg" und nichts hält es auf.