Runge-Kutta 4 implemetation blowing up
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Hello.
I am trying to solve the ODE system shown bellow using Runge-Kutta melhod, to compare with the solution given by the function ode23tb, that I successfully used to solve the exact same problem. I double checked the calculations and they seem to be correct. How can I avoid the extremely large numbers that I am getting in the Runge-Kutta coeficients?
fc = 10e3;
h = 1/(2000*fc);
t = h:h:0.01;
vD = 3;
vG = 10;
iL = 10/6.6;
s = double(...
[vD*ones(1,length(t));
vG*ones(1,length(t));
iL*ones(1,length(t))]);
for i = 1:length(t)-1
k1 = f(t(i),s(:,i));
k2 = f(t(i)+h/2,s(:,i)+h*k1/2);
k3 = f(t(i)+h/2,s(:,i)+h*k2/2);
k4 = f(t(i)+h,s(:,i)+h*k3);
s(:,i+1) = s(:,i+1) + 1/6*(k1 + 2*k2 + 2*k3 + k4);
end
plot(s(1,:));
function dydt = f(t,y)
VDD = 10;
Ld = 120e-6;
Cg = 50e-12;
Cb = 40e-12;
RS = 50;
RL = 6.6;
vG = y(1); vD = y(2);
iD = 10*(1/2.*(vG-3)+1/20.*log(2*cosh(10*(vG-3)))).*(1+0.003.*vD).*tanh(vD);
vs = 3+20*sin(2*pi*10e3*t);
dydt = [
1/Cg*((vs-y(1))/RS-iD-y(3)-y(2)/RL);
1/Cg*(vs-y(1))/RS-(Cb+Cg)/(Cb*Cg)*(iD+y(3)+y(2)/RL);
(y(2)-VDD)/Ld
];
end
2 Kommentare
David Wilson
am 19 Jan. 2021
Bearbeitet: David Wilson
am 19 Jan. 2021
It seems you have quite a stiff system to solve, and you are solving it for a very, very long time. You have two options:
(1) Use a decent stiff integrator such as ode23s, and let it choose the step size. Even then, you will need a small step size, and it seems no real point to go further than tfinal = 3e-4. The rest is just repetition.
(2) If you insist on using the brain-dead RK4, then you will need an extremely small step size to prevent instability. For example ode45 required h in the order of 1e-10. Of course then, you have no realistic idea of the accuracy, even if it is stable.
Antworten (1)
James Tursa
am 19 Jan. 2021
s(:,i+1) = s(:,i) + 1/6*(k1 + 2*k2 + 2*k3 + k4); % changed s(:,i+1) to s(:,i) on rhs
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!