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

1. Feb 23, 2013

### sinisterstarr

1. The problem statement, all variables and given/known data

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.

2. Relevant 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;

3. 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.

2. Feb 23, 2013

### Ray Vickson

Re: linear two-equation system, two variables to second derivative in

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.

3. Feb 23, 2013

### sinisterstarr

Re: linear two-equation system, two variables to second derivative in

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

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

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

4. Feb 23, 2013

### SteamKing

Staff Emeritus
Re: linear two-equation system, two variables to second derivative in

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''

5. Feb 24, 2013

### Ray Vickson

Re: linear two-equation system, two variables to second derivative in

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?

Last edited: Feb 24, 2013
6. Feb 24, 2013

### sinisterstarr

Re: linear two-equation system, two variables to second derivative in

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?

7. Feb 24, 2013

### Ray Vickson

Re: linear two-equation system, two variables to second derivative in

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.

8. Feb 25, 2013

### SteamKing

Staff Emeritus
Re: linear two-equation system, two variables to second derivative in

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.