Runge-Kutta 4 implemetation blowing up
1 view (last 30 days)
Show older comments
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 Comments
David Wilson
on 19 Jan 2021
Edited: David Wilson
on 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.
Answers (1)
James Tursa
on 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
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!