How to do Runge Kutta 4 with a second order ode?
12 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Andrew
am 28 Okt. 2014
Beantwortet: Diye
am 27 Jul. 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 Kommentare
Akzeptierte Antwort
James Tursa
am 28 Okt. 2014
Bearbeitet: James Tursa
am 28 Okt. 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)
Weitere Antworten (1)
Diye
am 27 Jul. 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 Kommentare
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!