Numerically Solving a System of Differential Equations Using a First-Order Taylor Series Approximation

I don't need specific code corrected for me (nor do I have any to show currently), just some guidance (and to see if what I need is even possible). I have a system of second-order differential equations (y'',y',y,z'',z',z) and all of their initial values which I have printed below. All non-bold letters are known constants (the functions are y'',y',y,z'',z',z; all are a function of time. At the end is the approximation to be used.

*the first instance of g in the first equation (at the start of the square brackets) should be a J

I am required to plot y and z and continue the calculations for this approximation process until y = pi, where the approximations then stop. I would then like to take the y value at this point (pi) and the z value at this point and then plug them into this equation (g is not a function)

Once this done, I can change the parameter e in those equations to e+dm for dm some very small increment, and then repeat the exact same process up to some value. At the end of this process I should have be able to have graphs for y and z against time and a graph for v_p against the e+dm values (mass).

My question is: What are some of the functions in MatLab I will need to know of to complete this task and how exactly would something like this be structured? I'm assuming something like this is actual doable in MatLab, at least.

[Ameer Hamza]

Have you tried ode45()? You will need to introduce 2 extra variables to convert it into a form compatible with ode45.

[Daniel]

I did try to make some use of ode45 but I was struggling to find a way to make it compatible with what I'm trying to do. Could I ask what you mean by 2 extra variables?

Thanks!

[Ameer Hamza]

ode45() requires that the equations should be represented in first order derivate form. For example see the 2nd example on the documentation page: https://www.mathworks.com/help/matlab/ref/ode45.html

[Daniel]

Alright. I made an attempt at something like that, which I'll show below. It doesn't seem to work, I'm not sure I understand too well how to apply that process to this? I do understand the idea needed but, just applying it is another thing.

The errors I get are: "Input argument 't' might be unused, although a later one is used. Consider it replacing by ~." (On the same line as "function") "The variable assigned to r might be unused" (same line as r =)

function dMdt = TREmodel(t,M) % declaring model
y = M(1); % y
z = M(2); % z
p = M(3); % y'
r = M(4); % z'
dpdt = ((g*sin(y)*IA+IB*(g*sin(z)*cos(y-z)-ls*(z.^2)*sin(y-z)-lla*(p^2)*cos(y-z)*sin(y-z))))/(ID+IC+(mp*(lla).^2*(sin(y-z)).^2)); % the second derivative theta
drdt = ((lla  .*  (((p^2)  .*  sin(y-z) - dpdt  .*  cos(y-z)) - g  .*  sin(z))))/ ls; % the second derivative phi
dMdt = [dpdt;drdt]; % output column
end 
M0 = [O;P]; % initial values theta phi
tRange = [0 20]; % time interval
[tSol, MSol] = ode45(@TREmodel,tRange,M0); % solving

[Ameer Hamza]

It seems that you are working in the correct direction but your model function should also return four derivatives. Try my answer below to see the if it works.

Answer 1

[Ameer Hamza]

You will need to introduce extra variables to convert the 2nd order equations to first order. I have assumed the following variables

y_dot = y1
z_dot = z1

From this I have

y_dotdot = y1_dot
z_dotdot = z1_dot

put y1_dot in place of y_dotdot and z1_dot in place of z_dotdot to get a first order system.

In the following code, i have assumed all the constants to be 1 and also the initial condition y(0) = 1, remaining initial conditions are zero.

Save the following function in your MATLAB path

function dydt = odefun(t, y) % y = [y z y1 z1]

[a,b,c,d,e,f,g,h,J] = deal(1);
dydt = zeros(4,1);
dydt(1) = y(3); % y_dot = y1
dydt(2) = y(4); % z_dot = z1
dydt(3) = (J*(a*b-c*d-e*f)*sin(y(1)) + e*f*J*sin(dydt(2))*sin(y(1)-y(2)) - f*dydt(1).^2*cos(y(1)-y(2))*sin(y(1)-y(2)))/...
  (a*b^2 + 1/12*c*(f+b)^2 + c*h^2 + e*f^2*sin(y(1)-y(2))^2);
dydt(4) = (f*(dydt(1)^2*sin(y(1)-y(2)) - dydt(3)*cos(y(1)-y(2))) - J*sin(y(2)))/...
  g;
end

and then run the following line;

[t,y] = ode45(@odefun, [0 pi], [1 0 0 0]);
plot(t,y,'-o')

The y will contain four columns, first for values of y and second for values of z, whereas 3rd and 4th are for the value of y1 and z1.

[Daniel]

Oh wow! Legend

I put all of this into my code but there's still for some reason an error. I still have an error relating to "Input argument 't' might be unused, although a later one is used. Consider it replacing by ~." (On the same line as "function") and I also have an error

"Error in trebuchetmate (line 26) y_dot = y1;"

Here is what I have put in. (The constants are now different from a, b, c, etc.. Sorry if that's somewhat confusing)

y_dot = y1;
z_dot = z1;
y_dotdot = y1_dot;
z_dotdot = z1_dot;
[t, y] = ode45(@odefun, [0 pi], [O P J K]);
plot(t,y,'-o');
function dydt = odefun(t, y) % y = [y z y1 z1]
dydt = zeroes(4,1);
dydt(1) = y(3); % y_dot = y1
dydt(2) = y(4); % z_dot = z1
dydt(3) = ((g*sin(y(1))*IA+IB*(g*sin(y(2))*cos(y(1)-y(2))-ls*(y^2)*sin(y(1)-y(2))-lla*(dydt(1).^2)*cos(y(1)-y(2))*sin(y(1)-y(2)))))/(ID+IC+(mp*(lla).^2*(sin(y(1)-y(2)))^2)); 
dydt(4) = ((lla  .*  (((dydt(1)^2)  .*  sin(y(1)-y(2)) - dydt(3)  .*  cos(y(1)-y(2))) - g  .*  sin(y(2)))))/ ls;
end

I can't really see what's different from yours? Perhaps I am misinterpreting some of your steps; sorry this process is very new to me.

[Ameer Hamza]

The first thing you mentioned is a warning, not error. You can suppress it by writing it like this

function dydt = odefun(~, y) % y = [y z y1 z1]

Also, you don't need to write

y_dot = y1;
z_dot = z1;
y_dotdot = y1_dot;
z_dotdot = z1_dot;

I wrote those lines just for explanation. Delete these lines and then run the code.

[Daniel]

Oh alright. Well thanks to you, I got that part to work! Now I just have to go learn all the other things required :x, which from what I've been looking at involves the event function.

Thanks so much! My hero. Here's the final code and graph if you're curious:

[t, y] = ode45(@odefun, [0 pi], [O P J K]);
plot(t,y,'-o');
legend('theta','phi','dtheta/dt','dphi,dt')
function dydt = odefun(~, y) % y = [y z y1 z1]
dydt = zeros(4,1);
g = 9.81; % gravity m/s^2
mcw = 16; % counterweight kg
lsa = 0.2; % length of short arm m
lla = 0.75; % length of long arm m
ls = 0.49; % length of sling m
h = 0.19 + 0.42; % height m
hf = 0.42; % height to fulcrum m
m = 0.42; % mass of arm kg
lcg = (lla - lsa)/2; % length to center of gravity
icw = 11.84 .*  mcw / 12; % inertia of counterweight kgm^2
mp = 0.05; % initial projectile mass kg
dm = 0.0005; % mass increments kg
dt = 0.0005; % step size s
O = acos((h-hf)/lsa); % initial theta value
P = 0; % initial phi value
OAZ = O - P; % the difference between O and P
J = 0; % initial angular velocity for theta value
K = 0; % initial angular velocity for phi value
IA = mcw * lsa - m * lcg - mp * lla; % making it easier to express constants
IB = mp * lla; % making it easier to express constants
IC = 1/12 * m * (lla + lsa)^2 + m * (lcg)^2; % making it easier to express constants
ID = mcw * (lsa)^2; % making it easier to express constants
U = ((g*sin(O)*IA+IB*(g*sin(P)*cos(OAZ)-ls*(K.^2)*sin(OAZ)-lla*(J.^2)*cos(OAZ)*sin(OAZ))))/(ID+IC+(mp*(lla).^2*(sin(OAZ)).^2)) ; %initial angular acceleration for theta value
N = ((lla  .*  (((J.^2)  .*  sin(OAZ) - U  .*  cos(OAZ)) - g  .*  sin(P))))/ ls; % initial angular acceleration for phi value
dydt(1) = y(3); % y_dot = y1
dydt(2) = y(4); % z_dot = z1
dydt(3) = ((g*sin(y(1))*IA+IB*(g*sin(y(2))*cos(y(1)-y(2))-ls*dydt(2).^2*sin(y(1)-y(2))-lla*dydt(1).^2*cos(y(1)-y(2))*sin(y(1)-y(2)))))/(ID + IC + mp*(lla).^2*(sin(y(1)-y(2)))^2);
dydt(4) = ((lla*((dydt(1).^2*sin(y(1)-y(2))-dydt(3)*cos(y(1)-y(2)))-g*sin(y(2)))))/ls;
end