Help with Difference Equation in MATLAB

Hi, all

I have a structural engineering book dated back at 1979 that presents a fairly odd equation to compute rotations at member joints.

This is:

Where C, V , h and K are constant properties. Theta is the joint rotation at the level considered, the level above and the level below (hence the subscript).

The book mentions a several ways of solving it. One being iteration procedures. It also refers to this equation as a difference equation, which I have no experience with at all. Although I have read that they are fairly similar to differential equations which i have some experience with.

My question is, is there a function or command in MATLAB that could be employed to solve this equation? Like there is for differential equations symbolically and numerically. And could somebody help me with this?

Many thanks,

Scott

[Scott Banks]

Thank you both for you assistance on this.

I am sorry to burden you further, but I have scanned two pages from the book. It might help if you have a read.

Torsten, if you notice on the first page under equation 17 it referes to equation 18 that it is a recurrence equation.

On the second page from equation 22 ownards I haven't mentioned anything because I am pretty clueless as to what it actually means. I don't know if you guys can make sense of it?

EDIT: C is the flexural stiffness of the columns and K is 6*G + C + C, where G is the flexural stiffness of the girders/beams. KGm is C/(6*G + C + C) and KGt = C/(6*G + C).

[Torsten]

(18) defines a linear system of equations with a tridiagonal coefficient matrix. This is confirmed by what is written in the first paragraph of 5.2 .

The paper is very old - thus solution of linear systems of equations was not standard at the time when it was published.

So if you write the equations in the form

A*theta = b

with A being a 8x8 coeffient matrix, b being a 8x1 vector and theta = [theta(1);...;theta(8)] being the solution vector and solve for theta as

theta = A\b

you will get your solution.

I did this for the equations you formulated. So if they are correct, you are done.

If you want to call (18) a recurrence relation, you can do this. The solution of this recurrence is done by solving a linear system of equations.

[Scott Banks]

Okay, thanks, Torsten. I will accept your answer!

Answer 1

[Torsten]

If V_i, h_i, C, K and theta_0 and theta_1 are given, you can solve for theta_(i+1) and compute theta in a loop.

theta(i+1) = ((V(i)*h(i)+V(i+1)*h(i+1))/4 - K*theta(i) + C*theta(i-1))/(-C)

[Scott Banks]

Hi Torsten,

I have tried that way before. The answers coverges very quickly. Look here:

% Structural proprties
E = 199947961.502;
Ic = 1.07E-03;
h = 3;

% Structural parameters
C = 1/2*4*E*Ic/h;
G = (2*E*Ic/5 + 1*E*Ic/3)*13;
K = (6*G + C);

thetaI = (30*h + 15*h)/(24*G) % initial approximation for storey 6

thetaI = 2.7579e-06

thetab = (45*h + 30*h)/(24*G) % initial approximation for storey 5

thetab = 4.5965e-06

for i = 1:10
   theta(i+1) = (15*h)/(4*C) + K/C*thetaI - 1*thetab  % Calculate theta for storey 7
end

theta = 1×2

1.0e-03 * 0 0.3137

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×3

1.0e-03 * 0 0.3137 0.3137

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×4

1.0e-03 * 0 0.3137 0.3137 0.3137

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×5

1.0e-03 * 0 0.3137 0.3137 0.3137 0.3137

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×6

1.0e-03 * 0 0.3137 0.3137 0.3137 0.3137 0.3137

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×7

1.0e-03 * 0 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×8

1.0e-03 * 0 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×9

1.0e-03 * 0 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×10

1.0e-03 * 0 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×11

1.0e-03 * 0 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137 0.3137

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

I randomly set the loop to 10 interations and as you can see to 4 decimal places convergence happens immediately.

[Torsten]

Your difference equation solution for theta seems to diverge:

% Structural proprties
E = 199947961.502;
Ic = 1.07E-03;
h = 3;
% Structural parameters
C = 1/2*4*E*Ic/h;
G = (2*E*Ic/5 + 1*E*Ic/3)*13;
K = (6*G + C);
thetaI = (30*h + 15*h)/(24*G) % initial approximation for storey 6

thetaI = 2.7579e-06

thetaI = 2.7579e-06

thetaI = 2.7579e-06

thetab = (45*h + 30*h)/(24*G) % initial approximation for storey 5

thetab = 4.5965e-06

thetab = 4.5965e-06

thetab = 4.5965e-06

theta(1) = thetaI;
theta(2) = thetab;
for i = 2:10
   theta(i+1) = (15*h)/(4*C) + K/C*theta(i) - 1*theta(i-1)  % Calculate theta for storey 7
end

theta = 1×3

1.0e-03 * 0.0028 0.0046 0.4751

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×4

0.0000 0.0000 0.0005 0.0413

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×5

0.0000 0.0000 0.0005 0.0413 3.5855

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×6

0.0000 0.0000 0.0005 0.0413 3.5855 311.1822

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×7

1.0e+04 * 0.0000 0.0000 0.0000 0.0000 0.0004 0.0311 2.7007

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×8

1.0e+06 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0003 0.0270 2.3439

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×9

1.0e+08 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0003 0.0234 2.0342

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×10

1.0e+10 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0002 0.0203 1.7655

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

theta = 1×11

1.0e+12 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0002 0.0177 1.5322

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Scott Banks]

Sorry, Torstem, I realised what I posted was rubbish.

This is how I have tried to solve the equation before. I condsider the entire frame and go joint by joint.

E = 2.0E+08;
I = 1.07E-03;
h = 3;
A = 0.0206;
C = 1/2*4*E*I/h;
G1 = E*(2*I/5 + I/3)*13;
Q = 15;

KGm = 6*G1 + C + C;
KGt = 6*G1 + C;

x1(1) = (Q*h/4)/(24*G1); % Storey 7 approximation
x2(1) = (2*Q*h/4 + Q*h/4)/(24*G1); % Storey 6 approximation
x3(1) = (3*Q*h/4 + 2*Q*h/4)/(24*G1); % Storey 5 approximation
x4(1) = (4*Q*h/4 + 3*Q*h/4)/(24*G1); % Storey 4 approximation
x5(1) = (5*Q*h/4 + 4*Q*h/4)/(24*G1); % Storey 3 approximation
x6(1) = (6*Q*h/4 + 5*Q*h/4)/(24*G1); % Storey 2 approximation
x7(1) = (7*Q*h/4 + 6*Q*h/4)/(24*G1 + 2*C); % Storey 1 approximation
x8(1) = 0;

% Start iterations....
for i = 1:10
     x1(i+1) = Q*h/(4*KGt) + C/KGm*x2(i);
     x2(i+1) = (Q*h + 2*Q*h)/(4*KGm) + C/KGt*x1(i) + C/KGm*x3(i);
     x3(i+1) = (2*Q*h + 3*Q*h)/(4*KGm) + C/KGm*x2(i) + C/KGm*x4(i);
     x4(i+1) = (3*Q*h + 4*Q*h)/(4*KGm) + C/KGm*x3(i) + C/KGm*x5(i);
     x5(i+1) = (4*Q*h + 5*Q*h)/(4*KGm) + C/KGm*x4(i) + C/KGm*x6(i);
     x6(i+1) = (5*Q*h + 6*Q*h)/(4*KGm) + C/KGm*x5(i) + C/KGm*x7(i);
     x7(i+1) = (6*Q*h + 7*Q*h)/(4*KGm) + C/KGm*x6(i);
     x8(i+1) = 0;
end

% create vector of rotations theta.
theta = [x1(end);x2(end);x3(end);x4(end);x5(end);x6(end);x7(end);x8(end)]

theta = 8×1

1.0e-04 * 0.0094 0.0276 0.0460 0.0643 0.0827 0.1011 0.1179 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

In this I made theta(i) the subject

Does this look any better you? have I done it correctly in your opinion?

[Torsten]

Your code has nothing to do with a recurrence relation. It seems that you try to solve a linear system of equations in the unknowns x1 - x8 by fixed point iteration.

You can solve it directly without iteration using the following code:

E = 2.0E+08;
I = 1.07E-03;
h = 3;
A = 0.0206;
C = 1/2*4*E*I/h;
G1 = E*(2*I/5 + I/3)*13;
Q = 15;
KGm = 6*G1 + C + C;
KGt = 6*G1 + C;

M =[1,-C/KGm,0,0,0,0,0,0;...
    -C/KGt,1,-C/KGm,0,0,0,0,0;...
    0,-C/KGm,1,-C/KGm,0,0,0,0;...
    0,0,-C/KGm,1,-C/KGm,0,0,0;...
    0,0,0,-C/KGm,1,-C/KGm,0,0;...
    0,0,0,0,-C/KGm,1,-C/KGm,0;...
    0,0,0,0,0,-C/KGm,1,0;...
    0,0,0,0,0,0,0,1];
rhs = [Q*h/(4*KGt);...
      (Q*h + 2*Q*h)/(4*KGm);...
      (2*Q*h + 3*Q*h)/(4*KGm);...
      (3*Q*h + 4*Q*h)/(4*KGm);...
      (4*Q*h + 5*Q*h)/(4*KGm);...
      (5*Q*h + 6*Q*h)/(4*KGm);...
      (6*Q*h + 7*Q*h)/(4*KGm);...
      0];
x = M\rhs

x = 8×1

1.0e-04 * 0.0094 0.0276 0.0460 0.0643 0.0827 0.1011 0.1179 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

I can't tell you if it's correct because I don't understand the underlying problem.

Answer 2

[John D'Errico]

Since V_i and h_i appear to vary with the index i, there is no analytical solution. (If V and h were constants, not changing with subscript, then an analytical solution will exist, and indeed, it would be not unlike the solution of a second order ODE.)

However, you can simply write a loop. You MUST know theta_0, and theta_1. Rewrite the equation, isolating theta_(i+1) on the left hand side. Now everything is known. Just compute the value of theta_(i+1) iteratively. This is no different from computing the Fibonacci numbers. Start at the beginning, and just write a loop.

If V and h are actually fixed, scalar constants that do NOT vary with the index, then the problem does indeed have a solution, and there are multiple ways you can solve for theta in an analytical form. But I won't go into the depth of explaining the solution for a case that I don't think exists at this time.

[Scott Banks]

Hi John, yes, they are fixed!

V is the shear at each joint at each storey, and h is simply the storey height. Here is a diagram that will help you visualise the problem:

'Q' here is what 'V' is and 'x 'here would be the 'theta'. The storey heights here is the same for each storey. You could assign 'h' here say 3 metres. Moreover, Q could be any value for the horizontal shear, lets say 15kN. Therefore, the total shear from top to bottom would be [15 30 45 60 75 90 105] kN.

[John D'Errico]

No. You don't understand. If V varies, across each level, then there is no analytical solution. Only if V has the same value for EACH level, can you find a simple analytical solution.

You have said that h is the same for each, but you have stated that V does vary. That it varies in a simple way, at least in your example, might possibly help.

Just use a loop!

[Scott Banks]

What about if we were to consider one level at each time. So then, for example, level 6, we can say:

Torsten, tried a loop earlier and it diverged.

[Torsten]

As said, the index i is not a recursion or iteration index, but an index for the solution variable x_i where i goes from 1 to 8. At least if we assume that your last code incorporates what you want to compute.

[John D'Errico]

This is a basic Gauss-Seidel iteration scheme. As you have wrritten it, if I recall properly, it will converge if the corresponding matrix probem is diagonally dominant. That is, you will need

1 > abs(C/KGM) + abs(C/KGt)

So what do you have?

E = 2.0E+08;
I = 1.07E-03;
h = 3;
A = 0.0206;
C = 1/2*4*E*I/h;
G1 = E*(2*I/5 + I/3)*13;
Q = 15;
KGm = 6*G1 + C + C;
KGt = 6*G1 + C;

Now, what are the ratios C/KGm and C/KGt?

C/KGm

ans = 0.0114

C/KGt

ans = 0.0115

Which means the iteration scheme you propose will be (fairly rapidly) convergent. As easily however, Torsten also showed how to solve the problem using backslash.