How to do Runge Kutta 4 with a second order ode?
12 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Risposta accettata
James Tursa
il 28 Ott 2014
Modificato: James Tursa
il 28 Ott 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)
Più risposte (1)
Diye
il 27 Lug 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 Commenti
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!