Help with Difference Equation in MATLAB
Mostrar comentarios más antiguos
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
3 comentarios
Scott Banks
el 17 de Ag. de 2025
Movida: Star Strider
el 17 de Ag. de 2025
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).
(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
el 19 de Ag. de 2025
Respuesta aceptada
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)
4 comentarios
Scott Banks
el 16 de Ag. de 2025
Editada: Scott Banks
el 16 de Ag. de 2025
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
thetab = (45*h + 30*h)/(24*G) % initial approximation for storey 5
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
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)]
In this I made theta(i) the subject
Does this look any better you? have I done it correctly in your opinion?
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
I can't tell you if it's correct because I don't understand the underlying problem.
Más respuestas (1)
John D'Errico
el 16 de Ag. de 2025
0 votos
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.
5 comentarios
Scott Banks
el 16 de Ag. de 2025
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
el 16 de Ag. de 2025
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
el 17 de Ag. de 2025
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.
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
C/KGt
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.
Categorías
Más información sobre Eigenvalue Problems en Centro de ayuda y File Exchange.
el 16 de Ag. de 2025
el 19 de Ag. de 2025
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
