MATLAB: Store alternating variable for for loop

alternating variablefor loop

I have this code that calulates Be(cidx) for each element by calculating B four times at 4 combinations of S1 and S2. Hence in the program B is calculated 16 times, 4 times for each element. I simply want to store B at each of the 16 times for each Be(cidx). I marked B with an annotation below in the code

E = 200*10^9;
v = 0.3;
t = 0.002;
f = 4000;
angle = 0;
%% D Matrix
D = (E/(1-v^2))*[1,v,0;v,1,0;0,0,0.5-0.5*v];
%% for loop
Ke1 = zeros (16,16);
Ke2 = zeros (16,16);
Ke3 = zeros (16,16);
Ke4 = zeros (16,16);
Be1 = zeros (3,16);
Be2 = zeros (3,16);
Be3 = zeros (3,16);
Be4 = zeros (3,16);
c1 = [0,0;0.02,0;0.02,0.01;0.01,0.01;0.01,0;0.02,0.005;0.015,0.01;0.005,0.005];
c2 = [0.02,0;0.03,0;0.03,0.01;0.02,0.01;0.025,0;0.03,0.005;0.025,0.01;0.02,0.005];
c3 = [0,0;0.01,0.01;0.01,0.025;0,0.025;0.005,0.005;0.01,0.0175;0.005,0.025;0,0.0125];
c4 = [0,0.025;0.01,0.025;0.01,0.04;0,0.04;0.005,0.025;0.01,0.0325;0.005,0.04;0,0.0325];
c = {c1, c2, c3, c4};
Ke = {Ke1, Ke2, Ke3, Ke4};
Be = {Be1, Be2, Be3, Be4};
for S1 = [-1/3^0.5,1/3^0.5];
    for S2 = [-1/3^0.5,1/3^0.5];
        dN1dS1 = (((1-S2)*(S2+2*S1))/4);
        dN1dS2 = (((1-S1)*(S1+2*S2))/4);
        dN2dS1 = (((1-S2)*(2*S1-S2))/4);
        dN2dS2 = (((1+S1)*(2*S2-S1))/4);
        dN3dS1 = (((1+S2)*(2*S1+S2))/4);
        dN3dS2 = (((1+S1)*(2*S2+S1))/4);
        dN4dS1 = (((1+S2)*(2*S1-S2))/4);
        dN4dS2 = (((1-S1)*(2*S2-S1))/4);
        dN5dS1 = ((-S1*(1-S2)));
        dN5dS2 = ((-1*(1-S1^2))/2);
        dN6dS1 = ((1-S2^2)/2);
        dN6dS2 = (-S2*(1+S1));
        dN7dS1 = (-S1*(1+S2));
        dN7dS2 = ((1-S1^2)/2);
        dN8dS1 = ((-1*(1-S2^2))/2);
        dN8dS2 = (-S2*(1-S1));
        H = 1;        
        N1 = [dN1dS1,dN2dS1,dN3dS1,dN4dS1,dN5dS1,dN6dS1,dN7dS1,dN8dS1;dN1dS2,dN2dS2,dN3dS2,dN4dS2,dN5dS2,dN6dS2,dN7dS2,dN8dS2]; 
        B1 = [1,0,0,0;0,0,0,1;0,1,1,0];
        B3 = [dN1dS1,0,dN2dS1,0,dN3dS1,0,dN4dS1,0,dN5dS1,0,dN6dS1,0,dN7dS1,0,dN8dS1,0;dN1dS2,0,dN2dS2,0,dN3dS2,0,dN4dS2,0,dN5dS2,0,dN6dS2,0,dN7dS2,0,dN8dS2,0;0,dN1dS1,0,dN2dS1,0,dN3dS1,0,dN4dS1,0,dN5dS1,0,dN6dS1,0,dN7dS1,0,dN8dS1;0,dN1dS2,0,dN2dS2,0,dN3dS2,0,dN4dS2,0,dN5dS2,0,dN6dS2,0,dN7dS2,0,dN8dS2];
        
        for cidx = 1 : 4 
               J = N1*c{cidx};
               Jdet = det(J);
               B2 = [J(2,2)/Jdet,-J(1,2)/Jdet,0,0;-J(2,1)/Jdet,J(1,1)/Jdet,0,0;0,0,J(2,2)/Jdet,-J(1,2)/Jdet;0,0,-J(2,1)/Jdet,J(1,1)/Jdet];
               B = B1*B2*B3; %% alternating variable - I want to store 16 times (four times for each Be (cidx))
               Be{cidx} = Be{cidx}+ B;
               k = H*t*B'*D*B*Jdet; 
               Ke{cidx} = Ke{cidx}+k;
        end
    end
end
Ke1 = Ke{1};
Ke2 = Ke{2};
Ke3 = Ke{3};
Ke4 = Ke{4};
Be1 = Be{1};
Be2 = Be{2};
Be3 = Be{3};
Be4 = Be{4};
%% Global Stifness Matrix Assembly
KGe1 = zeros(46,46);
KGe1 ([1,2,7,8,13,14,15,16,21,22,23,24,25,26,27,28],[1,2,7,8,13,14,15,16,21,22,23,24,25,26,27,28]) = Ke1 ([1,2,7,8,3,4,5,6,9,10,11,12,13,14,15,16],[1,2,7,8,3,4,5,6,9,10,11,12,13,14,15,16]);
KGe2 = zeros(46,46);
KGe2 ([3,4,5,6,13,14,15,16,23,24,29,30,31,32,33,34],[3,4,5,6,13,14,15,16,23,24,29,30,31,32,33,34]) = Ke2 ([3,4,5,6,1,2,7,8,15,16,9,10,11,12,13,14],[3,4,5,6,1,2,7,8,15,16,9,10,11,12,13,14]);
KGe3 = zeros(46,46);
KGe3 ([1,2,7,8,17,18,19,20,27,28,35,36,37,38,39,40],[1,2,7,8,17,18,19,20,27,28,35,36,37,38,39,40]) = Ke3 ([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16],[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]);
KGe4 = zeros(46,46);
KGe4 ([9,10,11,12,17,18,19,20,37,38,41,42,43,44,45,46],[9,10,11,12,17,18,19,20,37,38,41,42,43,44,45,46]) = Ke4 ([5,6,7,8,3,4,1,2,9,10,11,12,13,14,15,16],[5,6,7,8,3,4,1,2,9,10,11,12,13,14,15,16]);
KGlobal = KGe1+KGe2+KGe3+KGe4;
%% Gaussian Elimination
KBC = zeros(40,40);
KBC ([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40],[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40]) = KGlobal ([1,2,3,4,5,6,7,8,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,45,46],[1,2,3,4,5,6,7,8,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,45,46]);
F = [0;0;(f/6)*cosd(angle);(f/6)*sind(angle);(f/6)*cosd(angle);(f/6)*sind(angle);0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;(2*f/3)*cosd(angle);(2*f/3)*sind(angle);0;0;0;0;0;0;0;0;0;0;0;0];
U = inv(KBC)*F;
%% Post Processing
%% Element 1

%% Coordinates



x1 = 0; % node (1 - local / 1 - global)



y1 = 0; % node (1 - local / 1 - global)
x7 = 0.02; % node (2 - local / 7 - global)

y7 = 0; % node (2 - local / 7 - global)
x8 = 0.02; % node (3 - local / 8 - global)

y8 = 0.01; % node (3 - local / 8 - global)
x4 = 0.01; % node (4 - local / 4 - global)

y4 = 0.01; % node (4 - local / 4 - global)
x11= 0.01; % node (5 - local / 11 - global)

y11 = 0; % node (5 - local / 11 - global)
x12 = 0.02; % node (6 - local / 12 - global)

y12 = 0.005; % node (6 - local / 12 - global)
x13 = 0.015; % node (7 - local / 13 - global)

y13 = 0.01; % node (7 - local / 13 - global)
x14 = 0.005; % node (8 - local / 14 - global)

y14 = 0.005; % node (8 - local / 14 - global)
%% Element 2

%% Coordinates
x7 = 0.02; % node (1 - local / 7 - global)

y7 = 0; % node (1 - local / 7 - global)
x2 = 0.03; % node (2 - local / 2 - global)

y2 = 0; % node (2 - local / 2 - global)
x3 = 0.03; % node (3 - local / 3 - global)

y3 = 0.01; % node (3 - local / 3 - global)
x8 = 0.02; % node (4 - local / 8 - global)

y8 = 0.01; % node (4 - local / 8 - global)
x15 = 0.025; % node (5 - local / 15 - global)

y15 = 0; % node (5 - local / 15 - global)
x16 = 0.03; % node (6 - local / 16 - global)

y16 = 0.005; % node (6 - local / 16 - global)
x17 = 0.025; % node (7 - local / 17 - global)

y17 = 0.01; % node (7 - local / 17 - global)
x12 = 0.02; % node (8 - local / 12 - global)

y12 = 0.005; % node (8 - local / 12 - global)
%% Element 3

%% Coordinates
x1 = 0; % node (1 - local / 1 - global)
y1 = 0; % node (1 - local / 1 - global)
x4 = 0.01; % node (2 - local / 4 - global)

y4 = 0.01; % node (2 - local / 4 - global)
x9 = 0.01; % node (3 - local / 9 - global)

y9 = 0.025; % node (3 - local / 9 - global)
x10 = 0; % node (4 - local / 10 - global)

y10 = 0.025; % node (4 - local / 10 - global)
x14 = 0.005; % node (5 - local / 14 - global)

y14 = 0.005; % node (5 - local / 14 - global)
x18 = 0.01; % node (6 - local / 18 - global)

y18 = 0.0175; % node (6 - local / 18 - global)
x19 = 0.005; % node (7 - local / 19 - global)

y19 = 0.025; % node (7 - local / 19 - global)
x20 = 0; % node (8 - local / 20 - global)

y20 = 0.0125; % node (8 - local / 20 - global)
%% Element 4

%% Coordinates
x10 = 0; % node (1 - local / 10 - global)

y10 = 0.025; % node (1 - local / 10 - global)
x9 = 0.01; % node (2 - local / 9 - global)

y9 = 0.025; % node (2 - local / 9 - global)
x5 = 0.01; % node (3 - local / 5 - global)

y5 = 0.04; % node (3 - local / 5 - global)
x6 = 0; % node (4 - local / 6 - global)

y6 = 0.04; % node (4 - local / 6 - global)
x19 = 0.005; % node (5 - local / 19 - global)

y19 = 0.025; % node (5 - local / 19 - global)
x21 = 0.01; % node (6 - local / 21 - global)

y21 = 0.0325; % node (6 - local / 21 - global)
x22 = 0.005; % node (7 - local / 22 - global)

y22 = 0.04; % node (7 - local / 22 - global)
x23 = 0; % node (8 - local / 23 - global)

y23 = 0.0325; % node (8 - local / 23 - global)
%% Displacement for each element
%% Element 1
u1 = U(1,1); % NODE 1

v1 = U(2,1); % Node 1

u4 = U(7,1); % NODE 4

v4 = U(8,1); % Node 4

u7 = U(9,1); % NODE 7

v7 = U(10,1); % Node 7

u8 = U(11,1); % NODE 8

v8 = U(12,1); % Node 8

u11 = U(17,1); % NODE 11
v11 = U(18,1); % Node 11
u12 = U(19,1); % NODE 12

v12 = U(20,1); % Node 12

u13 = U(21,1); % NODE 13
v13 = U(22,1); % Node 13
u14 = U(23,1); % NODE 14

v14 = U(24,1); % Node 14

%% Strain & Stress



U_e1 = [u1;v1;u7;v7;u8;v8;u4;v4;u11;v11;u12;v12;u13;v13;u14;v14];
Strain_e1 = Be1*U_e1;
Stress_e1 = D*Strain_e1;
%% Element 2
u2 = U(3,1); % NODE 2
v2 = U(4,1); % Node 2
u3 = U(5,1); % NODE 3
v3 = U(6,1); % Node 3
u7 = U(9,1); % NODE 7
v7 = U(10,1); % Node 7
u8 = U(11,1); % NODE 8
v8 = U(12,1); % Node 8
u12 = U(19,1); % NODE 12
v12 = U(20,1); % Node 12
u15 = U(25,1); % NODE 15
v15 = U(26,1); % Node 15
u16 = U(27,1); % NODE 16
v16 = U(28,1); % Node 16
u17 = U(29,1); % NODE 17
v17 = U(30,1); % Node 17
%% Strain & Stress
U_e2 = [u7;v7;u2;v2;u3;v3;u8;v8;u15;v15;u16;v16;u17;v17;u12;v12];
Strain_e2 = Be2*U_e2;
Stress_e2 = D*Strain_e2;
%% Element 3
u1 = U(1,1); % NODE 1
v1 = U(2,1); % Node 1
u4 = U(7,1); % NODE 4
v4 = U(8,1); % Node 4
u9 = U(13,1); % NODE 9

v9 = U(14,1); % Node 9

u10 = U(15,1); % NODE 10

v10 = U(16,1); % Node 10

u14 = U(23,1); % NODE 14
v14 = U(24,1); % Node 14
u18 = U(31,1); % NODE 18
v18 = U(32,1); % Node 18
u19 = U(33,1); % NODE 19

v19 = U(34,1); % Node 19

u20 = U(35,1); % NODE 20
v20 = U(36,1); % Node 20
%% Strain & Stress
U_e3 = [u1;v1;u4;v4;u9;v9;u10;v10;u14;v14;u18;v18;u19;v19;u20;v20];
Strain_e3 = Be3*U_e3;
Stress_e3 = D*Strain_e3;
%% Element 4
u5 = 0; % NODE 5

v5 = 0; % NODE 5
u6 = 0; % NODE 6

v6 = 0; % NODE 6
u9 = U(13,1); % NODE 9
v9 = U(14,1); % Node 9
u10 = U(15,1); % NODE 10
v10 = U(16,1); % Node 10
u19 = U(33,1); % NODE 19
v19 = U(34,1); % Node 19
u21 = U(37,1); % NODE 21
v21 = U(38,1); % Node 21
u22 = 0; % NODE 22

v22 = 0; % NODE 22
u23 = U(39,1); % NODE 23
v23 = U(40,1); % Node 23
%% Strain & Stress
U_e4 = [u10;v10;u9;v9;u5;v5;u6;v6;u19;v19;u21;v21;u22;v22;u23;v23];
Strain_e4 = Be4*U_e4;
Stress_e4 = D*Strain_e4;

Best Answer

You asked this already, and I already showed you how to do it, with exact code using cell arrays, but you deleted that question.

B = {};
for S1 = [-1/3^0.5,1/3^0.5];
    for S2 = [-1/3^0.5,1/3^0.5];
        dN1dS1 = (((1-S2)*(S2+2*S1))/4);
        dN1dS2 = (((1-S1)*(S1+2*S2))/4);
        dN2dS1 = (((1-S2)*(2*S1-S2))/4);
        dN2dS2 = (((1+S1)*(2*S2-S1))/4);
        dN3dS1 = (((1+S2)*(2*S1+S2))/4);
        dN3dS2 = (((1+S1)*(2*S2+S1))/4);
        dN4dS1 = (((1+S2)*(2*S1-S2))/4);
        dN4dS2 = (((1-S1)*(2*S2-S1))/4);
        dN5dS1 = ((-S1*(1-S2)));
        dN5dS2 = ((-1*(1-S1^2))/2);
        dN6dS1 = ((1-S2^2)/2);
        dN6dS2 = (-S2*(1+S1));
        dN7dS1 = (-S1*(1+S2));
        dN7dS2 = ((1-S1^2)/2);
        dN8dS1 = ((-1*(1-S2^2))/2);
        dN8dS2 = (-S2*(1-S1));
        H = 1;        
        N1 = [dN1dS1,dN2dS1,dN3dS1,dN4dS1,dN5dS1,dN6dS1,dN7dS1,dN8dS1;dN1dS2,dN2dS2,dN3dS2,dN4dS2,dN5dS2,dN6dS2,dN7dS2,dN8dS2]; 
        B1 = [1,0,0,0;0,0,0,1;0,1,1,0];
        B3 = [dN1dS1,0,dN2dS1,0,dN3dS1,0,dN4dS1,0,dN5dS1,0,dN6dS1,0,dN7dS1,0,dN8dS1,0;dN1dS2,0,dN2dS2,0,dN3dS2,0,dN4dS2,0,dN5dS2,0,dN6dS2,0,dN7dS2,0,dN8dS2,0;0,dN1dS1,0,dN2dS1,0,dN3dS1,0,dN4dS1,0,dN5dS1,0,dN6dS1,0,dN7dS1,0,dN8dS1;0,dN1dS2,0,dN2dS2,0,dN3dS2,0,dN4dS2,0,dN5dS2,0,dN6dS2,0,dN7dS2,0,dN8dS2];
        
        for cidx = 1 : 4 
               J = N1*c{cidx};
               Jdet = det(J);
               Bidx = 
               B2 = [J(2,2)/Jdet,-J(1,2)/Jdet,0,0;-J(2,1)/Jdet,J(1,1)/Jdet,0,0;0,0,J(2,2)/Jdet,-J(1,2)/Jdet;0,0,-J(2,1)/Jdet,J(1,1)/Jdet];
               B{end+1} = B1*B2*B3; %% alternating variable - I want to store 16 times (four times for each Be (cidx))
               Be{cidx} = Be{cidx}+ B{end};
               k = H*t*B{end}'*D*B{end}*Jdet; 
               Ke{cidx} = Ke{cidx}+k;
        end
    end
end

Related Solutions

MATLAB: Problem with an if statement

Okay, think I know what you mean... this should do it.

c5 = [A(2:end,5); NaN];
idx = ((A(:,4)==1)|(A(:,4)==2))&(A(:,5)==c5)&(A(:,6)~=0);
B = A(idx,:);
A(idx,3) = abs(A(idx,6));

Tested with your example and gives the correct value for A.

MATLAB: How to solve for multiple variables using multiple equations

Here's a Newton-Raphson routine for the problem:

% Newton-Raphson Simultaneous Equations
% Six variables N1 to N6 need to be found from the six constraints f1 to f6.
% Variables N7 to N19 are found from N1 to N6.
% 
% The aim is to find values of the N that make each of the f's equal to
% zero (within some tolerance). 
% Set tolerances etc.
f = ones(6,1);  % large initial constraint values
tol = 10^-5;    % constraint tolerance
its = 0;        % iteration counter
flag = true;    % error testing flag
% Initial guesses for N1 to N6
N = [-0.5; 0.5; 0.005; 0.01; 0.2; 0.02];
% Newton-Raphson iteration loop
while flag && its < 60
    
    its = its+1;
    fold = f;
    Nold = N;
    N1 = N(1); N2 = N(2); N3 = N(3); N4 = N(4); N5 = N(5); N6 = N(6);
    k7 = 3.135;
    k8 = 0.41;
    k9 = 13.81;
    k10 = 76.1;
    k11 = 1.720E+05;
    k12 = 4.034E+04;
    k13 = 2.615E+04;
    k14 = 3.826E+04;
    k15 = 1.458E+03;
    k16 = 2.435E+04;
    k17 = 1.186E+08;
    k18 = 7.511E+02;
    k19 = 6.874E+02;
    N7 = k7*N1*N2;            dN7dN1 = k7*N2;         dN7dN2 = k7*N1; 
    N8 = k8*N2*N6;            dN8dN2 = k8*N6;         dN8dN6 = k8*N2; 
    N9 = k9*N1*N4;            dN9dN1 = k9*N4;         dN9dN4 = k9*N1;
    N10 = k10*N2^2*N3;        dN10dN2 = 2*k10*N2*N3;  dN10dN3 = k10*N2^2;
    N11 = k11*N2*N3^2*N5^3;   dN11dN2 = k11*N3^2*N5^3; dN11dN3 = 2*k11*N2*N3*N5^3; dN11dN5 = 3*k11*N2*N3^2*N5^2;
    N12 = k12*N3*N5^2;        dN12dN3 = k12*N5^2;     dN12dN5 = 2*k12*N3*N5;  
    N13 = k13*N2*N3^2*N5^2;   dN13dN2 = k13*N3^2*N5^2; dN13dN3 = 2*k13*N2*N3*N5^2; dN13dN5 = 2*k13*N2*N3^2*N5;
    N14 = k14*N1*N3*N5^2;     dN14dN1 = k14*N3*N5^2;  dN14dN3 = k14*N1*N5^2; dN14dN5 = 2*k14*N1*N3*N5;
    N15 = k15*N2*N3*N5;       dN15dN2 = k15*N3*N5;    dN15dN3 = k15*N2*N5;  dN15dN5 = k15*N2*N3;
    N16 = k16*N1*N3^2*N5;     dN16dN1 = k16*N3^2*N5;  dN16dN3 = 2*k16*N1*N3*N5;   dN16dN5 = k16*N1*N3^2;
    N17 = k17*N3^2*N5^3;      dN17dN3 = 2*k17*N3*N5^3; dN17dN5 = 3*k17*N3^2*N5^2;
    N18 = k18*N3*N5;          dN18dN3 = k18*N5;       dN18dN5 = k18*N3;
    N19 = k19*N3*N5;          dN19dN3 = k19*N5;       dN19dN5 = k19*N3;
    % Jacobian values
    % First row
    J11 = 1 + dN7dN1 + dN9dN1 + dN14dN1 + dN16dN1;
    J12 = 1 + dN7dN2 + dN8dN2 + dN10dN2 + dN11dN2 +dN13dN2 + dN15dN2;
    J13 = 1 + dN10dN3 + dN11dN3 + dN12dN3 + dN13dN3 + dN14dN3 + dN15dN3 + dN16dN3 + dN17dN3 + dN18dN3 + dN19dN3;
    J14 = 1 + dN9dN4;
    J15 = 1 + dN11dN5 + dN12dN5 + dN13dN5 + dN14dN5 + dN15dN5 + dN16dN5 + dN17dN5 + dN18dN5 + dN19dN5;
    J16 = 1 + dN8dN6;
    % Second row
    J21 = 0.1962*dN7dN1 -0.6575*(1 + dN7dN1 + dN9dN1 + dN14dN1 + dN16dN1);
    J22 = 0.1962*(1 + dN7dN2 + dN8dN2 + 2*dN10dN2 + dN11dN2 + dN13dN2 + dN15dN2) -0.6575*dN7dN2;
    J23 = 0.1962*(2*dN10dN3 + dN11dN3 +dN13dN3 + dN15dN3) - 0.6575*(dN14dN3 + dN16dN3);
    J24 = -0.6575*dN9dN4;
    J25 = 0.1962*(dN11dN5 + dN13dN5 + dN15dN5) - 0.6575*(dN14dN5 + dN16dN5);
    J26 = 0.1962*dN8dN6; 
    % Third row
    J31 = 0.0557*(dN14dN1 + 2*dN16dN1) - 0.4698*dN9dN1;
    J32 = 0.0557*(dN10dN2 + 2*dN11dN2 + 2*dN13dN2 + dN15dN2);
    J33 = 0.0557*(1 + dN10dN3 + 2*dN11dN3 + dN12dN3 + 2*dN13dN3 + dN14dN3 +dN15dN3 +2*dN16dN3 + 2*dN17dN3 + dN18dN3 + dN19dN3);
    J34 = -0.4698*(1+dN9dN4);
    J35 = 0.0557*(2*dN11dN5 + dN12dN5 + 2*dN13dN5 + dN14dN5 + dN15dN5 + 2*dN16dN5 + 2*dN17dN5 + dN18dN5 + dN19dN5);
    J36 = 0;
    % Fourth row
    J41 = 0.0099*(2*dN14dN1 + dN16dN1);
    J42 = 0.0099*(3*dN11dN2 + 2*dN13dN2 + dN15dN2) - 0.3525*dN8dN2;
    J43 = 0.0099*(3*dN11dN3 + 2*dN12dN3 +2*dN13dN3 + 2*dN14dN3 + dN15dN3 + dN16dN3 +3*dN17dN3 + dN18dN3 + dN19dN3);
    J44 = 0;
    J45 = 0.0099*(1 + 3*dN11dN5 + 2*dN12dN5 +2*dN14dN5 + dN16dN5 + 3*dN16dN5 + dN18dN5 + dN19dN5);
    J46 = -0.3525*(1 + dN8dN6);
    % Fifth row
    J51 = 0.0557*(1 + dN7dN1 + dN9dN1 + dN14dN1 + dN16dN1) - 0.1962*dN9dN1;
    J52 = 0.0557*dN7dN2;
    J53 = 0.0557*(dN14dN3 + dN16dN3);
    J54 = - 0.1962*(1 + dN9dN4);
    J55 = 0.0557*(dN14dN5 + dN16dN5);
    J56 = 0;
    % Sixth row
    J61 = 0.0099*dN7dN1;
    J62 = 0.0099*(1 + dN7dN2 + dN8dN2 + 2*dN10dN2 + dN11dN2 + dN13dN2 + dN15dN2) - 0.6575*dN8dN2;
    J63 = 0.0099*(2*dN10dN3 + dN11dN3 + dN13dN3 + dN15dN3);
    J64 = 0;
    J65 = 0.0099*(dN11dN5 + dN13dN5 + dN15dN5);
    J66 = 0.0099*dN8dN6 - 0.6575*(1 + dN8dN6);
    % Jacobian matrix
    J = [J11 J12 J13 J14 J15 J16;
         J21 J22 J23 J24 J25 J26;
         J31 J32 J33 J34 J35 J36;
         J41 J42 J43 J44 J45 J46;
         J51 J52 J53 J54 J55 J56;
         J61 J62 J63 J64 J65 J66];
    % Function values
    f1 = N1+N2+N3+N4+N5+N6+N7+N8+N9+N10+N11+N12+N13+N14+N15+N16+N17+N18+N19-1;
    f2 = 0.1962*(N2+N7+N8+2*N10+N11+N13+N15)-0.6575*(N1+N7+N9+N14+N16);
    f3 = .0557*(N3+N10+2*N11+N12+2*N13+N14+N15+2*N16+2*N17+N18+N19)-0.4698*(N4+N9);
    f4 = 0.0099*(N5+3*N11+2*N12+2*N13+2*N14+N15+N16+3*N17+N18+N19)-0.3525*(N6+N8);
    f5 = 0.0557*(N1+N7+N9+N14+N16)-0.1962*(N4+N9);
    f6 = 0.0099*(N2+N7+N8+2*N10+N11+N13+N15)-0.6575*(N6+N8);
    % Function vector
    f = [f1; f2; f3; f4; f5; f6];
    N = Nold - J\f;
    errmax = max(abs(f-fold));
    if errmax>tol
        flag = true;
    else
        flag = false;
    end
end
% Display resulting values of N1 to N6 and constraint values.
% N7 to N19 not displayed here, but appear in the workspace.
str = '       N        f';
disp(str)
disp([N f])

You should check it carefully!

Best Answer

Related Solutions

MATLAB: Problem with an if statement

MATLAB: How to solve for multiple variables using multiple equations

Related Question