Runge-Kutta 4 implemetation blowing up
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
David Wilson
el 19 de En. de 2021
Editada: David Wilson
el 19 de En. de 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.
Respuestas (1)
James Tursa
el 19 de En. de 2021
s(:,i+1) = s(:,i) + 1/6*(k1 + 2*k2 + 2*k3 + k4); % changed s(:,i+1) to s(:,i) on rhs
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!