How to do Runge Kutta 4 with a second order ode?
12 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Andrew
el 28 de Oct. de 2014
Respondida: Diye
el 27 de Jul. de 2018
I'm trying to do Runge Kutta with a second order ode, d2y/dx2+.5dy/dx+7y=0, with .5 step size from 0 to 5. My initial conditions are y'(0)=0 and y(0)=4.
This is my code so far:
s4=11; %Number of steps
h5=5/(s4-1); %Step size
x=linspace(0,5,s4);
y=[4; 0];
%Creating a for loop to define all the variables for Runge Kutta equation
for ii=1:length(x)-1
k1(:,ii)=ode1(x(ii),y(:,ii));
k2(:,ii)=ode1(x(ii)+.5*h5,y(:,ii)+k1(:,ii)*.5*h5);
k3(:,ii)=ode1(x(ii)+.5*h5,y(:,ii)+k2(:,ii)*.5*h5);
k4(:,ii)=ode1(x(ii)+h5,y(:,ii)+k3(:,ii)*h5);
y(2,ii+1)=y(2,ii)+1/6*(k1(2,ii)+2*k2(2,ii)+2*k3(2,ii)+k4(2,ii))*h5;
y(1,ii+1)=y(1,ii)+1/6*(k1(1,ii)+2*k2(1,ii)+2*k3(1,ii)+k4(1,ii));
end
I can't figure out why I'm not getting the right answer. Any ideas?
0 comentarios
Respuesta aceptada
James Tursa
el 28 de Oct. de 2014
Editada: James Tursa
el 28 de Oct. de 2014
What is the signature for ode1? I can't seem to find it in the doc. Where is your derivative function implemented? I don't see it in your code. If you look at the following link for RK4, the k's represent function derivatives at various predicted points, not integrated values as you seem to be doing.
E.g., I would have expected to see a derivative function something like this in your code somewhere:
function yp = y_deriv(x,y)
yp= [y(2);-7*y(1)-0.5*y(2)];
end
And this derivative function would be used to calculate the k's. But I don't see anything like this in your code. E.g.,
k1(:,ii)=y_deriv(x(ii),y(:,ii));
k2(:,ii)=y_deriv(x(ii)+.5*h5,y(:,ii)+k1(:,ii)*.5*h5);
k3(:,ii)=y_deriv(x(ii)+.5*h5,y(:,ii)+k2(:,ii)*.5*h5);
k4(:,ii)=y_deriv(x(ii)+h5,y(:,ii)+k3(:,ii)*h5);
Also looks like you are missing a *h5 on the last line. (you could combine the last two lines into one vector equation)
Más respuestas (1)
Diye
el 27 de Jul. de 2018
try solving it using this method.
firstly split the second order ode into systems of first order odes. let, dy/dx=z then, dz/dx=-5z-7y
in Matlab
function can %can is any name chosen for the script file
a0=[4 0]; [x, y]=ode45(@pol,[0:5:5],a0) y=a(:,1); z=a(:,2); plot(x,y) function dydx=pol(x,y) y=a(1); z=a(2); dydx=z; dydx=-5*z-7*y;
0 comentarios
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!