Linear two-equation system, two variables to second derivative in both

[sinisterstarr]

Homework Statement

I'm trying to solve a system of two second order linear differential equations with the ode45 function. It is a two degree of freedom problem with 2nd order derivatives of both variables, u and theta. I believe that's referred to as a "stiff matrix").

I'm very confident in the A and B matrices themselves, but not on the running of ode45. One of the two problems I'm having is the values in my y matrix. I'm just not sure what should be there. 1's have the same result as zeros. I know y is the dependent variable, I just don't understand what the function would want me to type there. What does its being zeros mean? What does its being ones mean?

My second problem is that I only know how to give a step input of force (f). As far as I can tell, the force input will be 628 Newtons the entire time. I tried to say f=628 when t<=0.015, but there is no "t" in the ode45 function. I'd like to give it an impulse if I can.

Homework Equations

My code so far is:

u0 = [0 ; 0];
theta0 = [0 ; 0];
FUN = @FUNA;
T = [0 10];
y0 = [u0 ; theta0];
[TT,YY] = ode45(FUN,T,y0);
plot(TT,YY(:,3));


function dy = FUNA(t,y)

%Constants of system
M = 70-5.876;
m = 5.876;
a = (((0.05)^2)+((0.13)^2))^0.5;
IG= 0.0233;
k = 500;
g = 9.81;
N = 2; %number of degrees of freedom

%Define sub matrices
MM=[(M+m) m*a ; m*a (IG+(m*a^2))];
K =[ 0 0 ; 0 (k-m*g*a)];
C =[0 0 ; 0 d*a*a];
F =[1; 0];
MI= inv(MM);
Z = zeros(N,N);
I = eye(N);

%Define matrices
A = [ Z I ; -MI*K -MI*C];
B = [ zeros(N,1); (MI*F) ];
f = 628;
y = [0; 0; 0; 0];

dy= A*y + B*f;

The Attempt at a Solution

I'd like to split this into four single order differential equations, but the second order derivative of both variables seems to stop me. I've run this code and I actually get plots, but they are either flat lines or ramp functions, depending on which of the four outputs I choose. I'm sure it's my use of the y matrix that's beating me.

 

[Ray Vickson]

Homework Statement

I'm trying to solve a system of two second order linear differential equations with the ode45 function. It is a two degree of freedom problem with 2nd order derivatives of both variables, u and theta. I believe that's referred to as a "stiff matrix").

I'm very confident in the A and B matrices themselves, but not on the running of ode45. One of the two problems I'm having is the values in my y matrix. I'm just not sure what should be there. 1's have the same result as zeros. I know y is the dependent variable, I just don't understand what the function would want me to type there. What does its being zeros mean? What does its being ones mean?

My second problem is that I only know how to give a step input of force (f). As far as I can tell, the force input will be 628 Newtons the entire time. I tried to say f=628 when t<=0.015, but there is no "t" in the ode45 function. I'd like to give it an impulse if I can.

Homework Equations

My code so far is:

u0 = [0 ; 0];
theta0 = [0 ; 0];
FUN = @FUNA;
T = [0 10];
y0 = [u0 ; theta0];
[TT,YY] = ode45(FUN,T,y0);
plot(TT,YY(:,3));


function dy = FUNA(t,y)

%Constants of system
M = 70-5.876;
m = 5.876;
a = (((0.05)^2)+((0.13)^2))^0.5;
IG= 0.0233;
k = 500;
g = 9.81;
N = 2; %number of degrees of freedom

%Define sub matrices
MM=[(M+m) m*a ; m*a (IG+(m*a^2))];
K =[ 0 0 ; 0 (k-m*g*a)];
C =[0 0 ; 0 d*a*a];
F =[1; 0];
MI= inv(MM);
Z = zeros(N,N);
I = eye(N);

%Define matrices
A = [ Z I ; -MI*K -MI*C];
B = [ zeros(N,1); (MI*F) ];
f = 628;
y = [0; 0; 0; 0];

dy= A*y + B*f;

The Attempt at a Solution

I'd like to split this into four single order differential equations, but the second order derivative of both variables seems to stop me. I've run this code and I actually get plots, but they are either flat lines or ramp functions, depending on which of the four outputs I choose. I'm sure it's my use of the y matrix that's beating me.

Why don't you first write out in detail the actual two DEs in the two variables? That way, we can understand what your problem is, without trying to decipher your post. Once we know the system we might be able to make useful suggestions.

 

[sinisterstarr]

Sorry, I guess that would be a lot easier! These are the two equations:

(m*a)u” + (I + m*a² )θ” + (d*a² )θ’ + (K - m*g*a)θ = 0

(M + m)u” + (m*a)θ” = -f₄

 

[SteamKing]

The numerical methods used to solve ODEs work only on first order equations. If you have second or higher order ODEs, then these equations must be converted into a system of first order equations.

For you equations, consider:

w = u'
w' = u''
alpha = theta'
alpha' = theta''

 

[Ray Vickson]

Sorry, I guess that would be a lot easier! These are the two equations:

(m*a)u” + (I + m*a² )θ” + (d*a² )θ’ + (K - m*g*a)θ = 0

(M + m)u” + (m*a)θ” = -f₄

This system is solvable symbolically (assuming that m, a, I, d, K, g, M, f4 are all constants), so there is no *need* to do it numerically.
Even for f4 = f4(t) = general function of t, you can write down the solution in closed-form (although it is messy). Is there a good reason why you want to do it numerically?

 

[sinisterstarr]

Thank you all for the advice!

SteamKing, would using this substitution mean that the non-differentiated theta variable in the first equation would become an integral of alpha? If I'm correct about that, I feel like I'd be in a worse spot. What should I do with that integral?

Ray Vickson, I'm not intent on solving the problem numerically. I'm actually trying to do it both ways. I've used the dsolve function and gotten two very long, messy symbolic solutions. Using those equations, if I put a single value of t (t=1), I'll get values for u and theta. However, when I try to plot the output with t = 0:0.01:10, I get error. Here's the attempt to solve it symbolically I've made:

M =70-5.876;
m =5.876;
a =(((0.05)^2)+((0.13^2))^0.5);
IG = 0.0233;
d = 500;
k = 500;
g = 9.81;
f = 628;

%y is u, x is theta

syms M m a IG d k g y(t) x(t)

eqn1 = (M+m)*diff(y,2) + M*diff(x,2) == -f;
eqn2 = m*a*diff(y,2) + (IG + m*a*a)*diff(x,2) + (d*a*a)*diff(x) + (k - m*g*a)*x == 0;

z = dsolve(eqn1,eqn2, y(0)==0, Dy(0)==0, x(0)==0, Dx(0)==0, 't');

z.x
z.y

Can you tell me please if I'm on the right path to solving it symbolically?

 

[Ray Vickson]

Thank you all for the advice!

SteamKing, would using this substitution mean that the non-differentiated theta variable in the first equation would become an integral of alpha? If I'm correct about that, I feel like I'd be in a worse spot. What should I do with that integral?

Ray Vickson, I'm not intent on solving the problem numerically. I'm actually trying to do it both ways. I've used the dsolve function and gotten two very long, messy symbolic solutions. Using those equations, if I put a single value of t (t=1), I'll get values for u and theta. However, when I try to plot the output with t = 0:0.01:10, I get error. Here's the attempt to solve it symbolically I've made:

M =70-5.876;
m =5.876;
a =(((0.05)^2)+((0.13^2))^0.5);
IG = 0.0233;
d = 500;
k = 500;
g = 9.81;
f = 628;

%y is u, x is theta

syms M m a IG d k g y(t) x(t)

eqn1 = (M+m)*diff(y,2) + M*diff(x,2) == -f;
eqn2 = m*a*diff(y,2) + (IG + m*a*a)*diff(x,2) + (d*a*a)*diff(x) + (k - m*g*a)*x == 0;

z = dsolve(eqn1,eqn2, y(0)==0, Dy(0)==0, x(0)==0, Dx(0)==0, 't');

z.x
z.y

Can you tell me please if I'm on the right path to solving it symbolically?

I don't know what system you are using, but in Maple 11 here is what I have (using symbol 'dd' instead of 'd', because 'd' is a reserved symbol in Maple):

params:={M =70-5.876,m =5.876,a =sqrt((0.05)^2+(0.13)^2),IG = 0.0233,
> dd = 500,k = 500,g = 9.81,f = 628 };
params := {M = 64.124, m = 5.876, a = 0.1392838828, IG = 0.0233,

dd = 500, k = 500, g = 9.81, f = 628}

eqn1:=(M+m)*diff(y(t),t$2) + M*diff(x(t),t$2) = -f: lprint(eqn1);(M+m)*diff(diff(y(t),t),t)+M*diff(diff(x(t),t),t) = -f
eqn2:= m*a*diff(y(t),t$2) + (IG + m*a*a)*diff(x(t),t$2) + (dd*a*a)*diff(x(t),t)+ (k - m*g*a)*x(t) =0: lprint(eqn2);m*a*diff(diff(y(t),t),t)+(IG+m*a^2)*diff(diff(x(t),t),t)+dd*a^2*diff(x(t),t)+(k-m*g*a)*x(t) = 0

sol:=dsolve({eqn1,eqn2,y(0)=0,D(y)(0)=0,x(0)=0,D(x)(0)=0},{x(t),y(t)}):


Note: the final result here is the symbolic solution for any unspecified values of the constants. It is lengthy and complicated, so there is little point in writing it out here. (It takes 3 pages of Maple output to print it.) It would have been much shorter if we had used the numerical values of the inputs before solving, but I just wanted to illustrate the fact that it could be done completely symbolically. Substituting the numerical values into the general formula gives:

sol1:=subs(params,sol): lprint(sol1);
{y(t) = Int(1.000000000*Int(4.013583972*exp(-21.50895008*_z1)+6.969036900*exp(37.34733541*_z1),_z1 = 0 .. _z1)-8.971428574*_z1,_z1 = 0 .. t), x(t) = -.1293011615e-6*(12237.52123*4^(1/2)+48768.43222)*exp(.1166302307e-1*(-1261.600126*4^(1/2)+679.0000004)*t)+.2245136944e-6*(12237.52123*4^(1/2)-48768.43222)*exp(.1166302307e-1*(1261.600126*4^(1/2)+679.0000004)*t)+.1492466504e-1}

yt:=rhs(sol1[1]): xt:=rhs(sol1[2]):
value(yt) assuming real: yt:=%:

lprint(xt);-.1293011615e-6*(12237.52123*4^(1/2)+48768.43222)*exp(.1166302307e-1*(-1261.600126*4^(1/2)+679.0000004)*t)+.2245136944e-6*(12237.52123*4^(1/2)-48768.43222)*exp(.1166302307e-1*(1261.600126*4^(1/2)+679.0000004)*t)+.1492466504e-1

lprint(yt);-.1367184601e-1+.1332761138e-10*t+.8675488339e-2*exp(-21.50895008*t)+.4996357667e-2*exp(37.34733541*t)-4.485714287*t^2

Now x(t) = xt and y(t) = yt plot perfectly well. They would also have plotted if I had not bothered to simplify y(t) using 'value(..) assuming real', but I wanted to have an explicit formula before plotting.

 

[SteamKing]

Any variables which are not differentiated, like theta, would remain as original in the equation. What you are trying to do with the substitutions is to reduce the degree of the differential equation to the first order. The values of variables like alpha' and w' are used only in the numerical calculations of the original ODEs. In other words, you treat alpha' and w' like 'dummy variables' in order to obtain a numerical solution.