I want help in writing a matlab code so that I wouldn't have to manually input the various steps as I have shown below. Thanks

E = 2900; % elastic modulus
A = 15;   % area
Ix = 19;  % modulus of inertia about x-axis
Iy = 11;  % modulus of inertia about y-axis
Ixy = 12; % modulus of inertia about x-y axis
K = [ E*A 1 1 1; 1 E*Ix 1 1; 1 1 E*Iy 1; 1 1 1 E*Ixy ]; % stiffness matrix
u = [ 4; 6; 7; 8 ]; % deformation vector
F = [ 15;29;34;45]; % force vector
Q = [ 100;78;85;79]; % maximum force vector
R = Q - F; % initialize R
dF = 2*ones(4,1); % incremental force vector
    T  = F + dF;  % new force vector    
    R  = Q - T;   % residual force vector
    m  = pinv(K)*R; % current deformation vector
    q  = u + m;       %% end of first step 
    T1 = T + dF;  % new force vector (2)
    R1 = Q - T1;   % residual force vector (2)
    K1 = T1\q;     % updated stiffness matrix (2)  
    m1 = pinv(K1)*R1; % updated current deformation vector (2)
    q1 = q + m1;      %% end of second step
    T2 = T1 + dF;  % new force vector (3)
    R2 = Q - T2;   % residual force vector (3)
    K2 = T2\q1;     % updated stiffness matrix (3)
    m2 = pinv(K2)*R2; % updated current deformation vector (3)
    q2 = q1 + m2;      %% end of third step
    T3 = T2 + dF;  % new force vector (4)
    R3 = Q - T3;   %residual force vector (4)
    K3 = T3\q2;    % updated stiffness matrix (4)
    m3 = pinv(K3)*R3; % updated current deformation vector (4)
    q3 = q2 + m3;     %% end of fourth step
    T4 = T3 + dF; % new force vector (5)
    R4 = Q - T4;  % residual force vector (5)
    K4 = T4\q3;   % updated stiffness matrix (5)
    m4 = pinv(K4)*R4; % updated current deformation vector (5)
    q4 = q3 + m4;     %% end of fifth step

Answer 1

[Walter Roberson]

Unfortunately, someone who was not familiar with the field would have little idea what is going on. What is K -- is that spring constants? Are Ix an Iy and Ixy momentums ? Is Q "quality factor" ?

You are not updating K during your loop, so it is not clear why you are calculating pinv(K) inside the loop instead of before the loop?

In order to "write clearly" then you need comments -- even if you were using variable names that explained more, you would at least have comments at the beginning explaining the purpose of the code.

[George Adomako Kumi]

I have edited it with the required comments. My original thoughtt is for everything to be calculated inside the loop.

[Walter Roberson]

    i = 0;      % iteration begins

That should not be inside the loop.

[Walter Roberson]

No comments about what T, R, m, q, S mean ?

[George Adomako Kumi]

i HAVE DULY COMMENTED IT AS REQUIRED

[Walter Roberson]

Commented version does not seem to be present ?

[George Adomako Kumi]

I have just commented it again.

Thanks

Answer 2

[John D'Errico]

R is a matrix. So the statemen

while R

is meaningless. As well, this if statement

if R <= 0.0001;

break

end

is also useless.

The problem is, there is no indication of what you are trying to do. What problem you are trying to solve. (I recall you have asked essentially the same question before.) So it is difficult to know if your implementation is correct.

If your goal is to iterate until the MATRIX R is essentially zero, then a good idea is to test the matrix norm, applied to R.

Next, you can put that test into the while statement. So there is absolutely no need to have an if statement ainside the loop, since the loop will iterate until that test in the while loop is not satisfied.

What else? It is often a bad idea to use i as a loop counter, since i is defined in MATLAB as sqrt(-1).

What else? This is illegal syntax:

dF = 2*ones*(4,1);

What else? R was never defined initially in your code, you you use it in the call to the while loop before it is defined inside the loopitself. You NEED to initialize R.

Given all of that, I might try this code:

E = 2900;
A = 15;
Ix = 19;
Iy = 11;
Ixy = 12;
K = [ E*A 1 1 1; 1 E*Ix 1 1; 1 1 E*Iy 1; 1 1 1 E*Ixy ];
u = [ 4; 6; 7; 8 ];
F = [ 15;29;34;45];
Q = [ 100;78;85;79];
R = Q - F; % initialize R
dF = 2*ones(4,1);
tol = 1e-8;
iter = 0;
while norm(R) > tol 
    iter = iter + 1
    T  = F + dF;
    R  = Q - T;
    m  = pinv(K)*R;
    q  = u + m;
    S  = F/q;
end

However, this code seems never to converge.

[George Adomako Kumi]

Thanks