RK4/AB4, need help with correct code for 2 second order equations in Matlab

6 views (last 30 days)
I have examples of code to follow for one second order equation but struggling to write correct code for 2 second order equations. grateful for any help!
  4 Comments
James Tursa
James Tursa on 4 Apr 2020
Let us know if you need more help. We can't give much advice until we see the differential equations and the code you have written.
Anne Ellaway
Anne Ellaway on 4 Apr 2020
Edited: darova on 6 Apr 2020
here's the code I have so far
function tryout()
close all
h=0.01;
t0=0;
tf=10;
N=(tf-t0)/h;
t=(t0:h:tf);
%create array
w=zeros(4,N+1); %x1
x=zeros(4,N+1); %x2
y=zeros(4,N+1); %v1
z=zeros(4,N+1); %v2
%initial conditions
w0(1)=1; %x1(0)=1
x0(2)=6; %x2(0)=6
y0(1)=0; %dx1/dt=v1=0
z0(2)=3; %dx2/dt=v2=3
%inputting functions
f_x1='f_x1';
f_x2='f_x2';
f_V1='f_V2';
f_V2='f_V2';
%runge kutta
for i=1:N
k1 = h * feval(f_x1, w(i), x(i), y(i), z(i));
k2 = h * feval(f_x1, w(i)+(k1/2), x(i)+(k1/2), y(i)+(k1/2), z(i)+(k1/2));
k3 = h * feval(f_x1, w(i)+(k2/2), x(i)+(k2/2), y(i)+(k2/2), z(i)+(k2/2));
k4 = h * feval(f_x1, w(i)+k3, x(i)+k3, y(i)+ k3, (z(i)+k3));
w(i+1) = w(i) + ( k1 + 2.*k2 + 2.*k3 + k4 ) ./ 6;
x(i+1) = x(i) + ( k1 + 2.*k2 + 2.*k3 + k4 ) ./ 6;
y(i+1) = y(i) + ( k1 + 2.*k2 + 2.*k3 + k4 ) ./ 6;
z(i+1) = z(i) + ( k1 + 2.*k2 + 2.*k3 + k4 ) ./ 6;

Sign in to comment.

Accepted Answer

James Tursa
James Tursa on 4 Apr 2020
Edited: James Tursa on 4 Apr 2020
So, first define a 4-element state vector. To keep the nomenclature the same as the MATLAB docs, I will use the variable name y. The elements of y are defined as
y(1) = x1
y(2) = x2
y(3) = x1dot
y(4) = x2dot
Next step is to solve your matrix differential equations for x1dotdot and x2dotdot so that you have two expressions. You can do this easily on paper. I.e., do the matrix multiply on paper and then solve the two equations for the highest order derivatives. You will get:
x1dotdot = some expression involving x1, x2, x1dot, and x2dot
x2dotdot = some expression involving x1, x2, x1dot, and x2dot
Then rewrite these expressions using the definitions above, so that you will get
x1dotdot = some expression involving y(1), y(2), y(3), and y(4)
x2dotdot = some expression involving y(1), y(2), y(3), and y(4)
Now you are ready to write code. The derivative function will be
dy = @(t,y) [y(3); y(4); the expression for x1dotdot; the expression for x2dotdot]
This can be used in your RK4 code, or even passed into ode45( ).
  10 Comments
James Tursa
James Tursa on 5 Apr 2020
The backslash operator was a bit of overkill here since the equations can be easily solved by hand. But basically you have this matrix equation
M*Xdotdot + K*X + C*Xdot = 0
where M, K, and C are as above and you have the following state vector definitions
X = [x1;x2]
Xdot = [x1dot;x2dot]
Xdotdot = [x1dotdot;x2dotdot]
You are trying to solve for the highest order derivatives in the equation, which is Xdotdot. So just doing the math:
M*Xdotdot = -K*X - C*Xdot
Now you have an equation in the form Ax=b. MATLAB solves this with the backslash operator, x = A\b. So in your case we get
Xdotdot = M \ (-K*X - C*Xdot)
Put in terms of your variable names, X = [x1;x2] = y(1:2), and Xdot = [x1dot;x2dot] = y(3:4).
Anne Ellaway
Anne Ellaway on 5 Apr 2020
Thank you so much again James, I really appreciate all your help in this. All the best to you in these challenging times.

Sign in to comment.

More Answers (0)

Categories

Find more on Function Creation 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!