I want to shorten this code but dont dont know how???

Question

Kundan 2023-9-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2023142-i-want-to-shorten-this-code-but-dont-dont-know-how

回答： Shreshth 2024-1-5

function [theta, M] = assignment(fc,fy,As,xl,xu,yl,yu,nx,tol,thetaf,m)
 % Input data, data can be either for under reinforced or over reinforced
 function [fc,fy,As,xl,xu,yl,yu,nx,tol,thetaf,m] = data
  fc = 30;             % Mean compressive strength of concrete
  fy = 415;            % Yield strength of steel
  As = [324 50 50; 324 150 50]; % [(Cross sectional area of rebar) (x-coordinate of rebar) (y-coordinate of rebar)] 
  xl = 0;              % Lower limit of x-coordinate
  xu = 200;            % Upper limit of x-coordinate
  yl = 0;              % Lower limit of y-coordinate
  yu = 500;            % Upper limit of y-coordinate
  nx = 200;            % No of grid points in each direction 
  tol = 1e-3;          % Tolerance for finding neutral axis
  thetaf = 0.08/1000;  % Final value of curvature
  m = 100;             % Number of increments
 end
 
 
 % Stress strain curve for concrete
 
 function sigma= sigmac (e, fc)
  ec2 = 2.0 / 1000 ;
  ecu = 3.5 / 1000 ;
  n = 2 ;
  if e < 0
   sigma = 0 ;
  else
   if e > ecu
    sigma = 0 ;
   elseif e > ec2
    sigma = fc ;
   else
    sigma = fc * (1 - (1 - e / ec2)^n) ;
   end
  end
 end
 % Stress strain curve for steel
 function sigma = sigmas (e, fy)
  Es = 200000;
  esy = fy / Es ;
  esu = 10 / 1000 ;
  if abs (e) > esu
   sigma = 0 ;
  elseif abs (e) > esy
   sigma = fy * e / abs (e) ;
  else
   sigma = e * Es ;
  end
 end
 % Finding the geometry of cross-section
 function I = geometry(x)
  xv = [0 200 200 0];
  yv = [0 0 500 500];
  I = inpolygon(x(1),x(2),xv,yv);
 end
 % Moment curvature diagram
 % {Evaluating curvature considering the cross-section and eliminating other
 % point which are outside the cross-section} 
 function [theta,M] = diagram(fc,fy,As,xl,xu,yl,yu,nx,tol,thetaf,m)
  A0 = (xu-xl)*(yu-yl);
  n = nx^2;
  w = A0/n;
  nr = size(As,1);     % Number of rebars
  du = 1/nx;
  u = (du / 2): du: (1-du / 2) ;
  k = 0 ;
  for i =1: nx
   for j = 1: nx
    k = k + 1 ;
    S(k,:) = [ u(i) u(j) ] ;
   end
  end
  S (:, 1) = xl + S (:, 1) * (xu - xl) ;
  S (:, 2) = yl + S (:, 2) * (yu - yl) ;
  for i =1:n
   I (i) = geometry (S(i,:)) ;
  end
  S(I ==0,:) = [ ] ;
  ns = size (S, 1) ;
  dtheta = thetaf / m;
  for i =2: m
   theta (i) = i * dtheta ;
   M(i) = momrot (theta(i), ns, S, w, nr, As, fc, fy, yl, yu, tol) ;
  end
 end 
 % Finding the position of neutral axis
 function M = momrot (theta, ns, S, w, nr, As, fc, fy, yl, yu, tol)
  a = yl ;
  b = yu ;
  Na = resultants (theta, ns, S, w, nr, As, a, fc, fy) ;
  h0 = b - a;
  it = ceil (log2 (h0 / tol)) ;
  for i = 1: it
   c = 0.5 * (a + b) ;
   Nc = resultants (theta, ns, S, w, nr, As, c, fc, fy) ;
   if Na * Nc < 0
    b = c ;
   else
    a = c ;
    Na = Nc ;
   end
  end
  yc = 0.5 * (a + b) ;
  [ N, M ] = resultants (theta, ns, S, w, nr, As, yc, fc, fy) ;
 end
 % Finding the moment
 function [ N, M ] = resultants (theta, ns, S, w, nr, As, yc, fc, fy)
  N = 0 ;
  M = 0 ;
  for i =1: ns
   y = S(i, 2) ;
   e = tan(theta) * (y - yc) ;
   sigma = sigmac (e, fc) ;
   N = N + w * sigma ;
   M = M + w * sigma * (y - yc) ;
  end
  for i =1: nr
   y = As (i, 3) ;
   e = tan(theta) * (y - yc) ;
   sigma = sigmas (e, fy) ;
   N = N + sigma * As (i, 1) ;
   M = M + sigma * As (i, 1) * (y - yc) ;
  end
 end
 [ fc, fy, As, xl, xu, yl, yu, nx, tol, thetaf, m ]=data ;
 [ theta, M ] = diagram(fc, fy, As, xl, xu, yl, yu, nx, tol, thetaf, m) ;
 max(M);
 disp(M);
 plot (theta, M);
 xlabel("Curvature (1/mm)");
 ylabel("Moment (Nmm)");
end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Mathieu NOE 2023-9-20

why not simply store individually your many functions into a toolbox (must be declared on your path)

then your main code needs only to call (in one line) your functions

no need to have them permanently present at the end of your script

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Shreshth 2024-1-5

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2023142-i-want-to-shorten-this-code-but-dont-dont-know-how#answer_1383501

在 MATLAB Online 中打开

Hello Kundan,

I understand that you want to shorten the given code. This can be achieved through several ways. Some of them I have mentioned below.

1.Vectorization: MATLAB is highly optimized for vector and matrix operations. Wherever possible, loops should be replaced with vectorized operations. This can significantly reduce the length of the code and often improve performance.

Original version:

for i =1: ns
  y = S(i, 2);
  e = tan(theta) * (y - yc);
  sigma = sigmac(e, fc);
  N = N + w * sigma;
  M = M + w * sigma * (y - yc);
End

Vectorized version:

e = tan(theta) * (S(:, 2) - yc);
sigma = arrayfun(@(strain) sigmac(strain, fc), e);
N = sum(w * sigma);
M = sum(w * sigma .* (S(:, 2) - yc));

2.Inline Functions: MATLAB allows the creation of anonymous functions (lambdas) for simple operations. These can be used to replace standalone functions that are only used once or have a very simple body.

Original standalone function:

function sigma = sigmac(e, fc)
  % ... function body ...
end

Inline version:

sigmac = @(e, fc) (e >= 0 & e <= 3.5/1000) .* fc * (1 - (1 - e / (2.0/1000)).^2);

3.Removing Redundant Code: If there are checks or initializations that are unnecessary or can be simplified, they should be removed or combined.

Original code checking for stress in concrete:

if e < 0
  sigma = 0;
else
  % ... more conditions ...
end

Simplified version (assuming e is always non-negative in context)

sigma = fc * (1 – (1 – e / ec2).^n) .* (e <= ec2) + (e > ec2 & e <= ecu) * fc;

4.Combining Functions: If functions are only called once and are specific to the problem at hand, their functionality can sometimes be combined into a single function to reduce the overhead of multiple function calls.

5.Use Built-in Functions: MATLAB has a wide range of built-in functions that can often replace custom code. For example, the ‘inpolygon’ function can be used to test if points are inside a polygon, replacing custom geometry functions.

Original custom geometry check:

function I = geometry(x)
  % ... function body ...
End

Using inpolygon:

I = inpolygon(S(:,1), S(:,2), [0, 200, 200, 0], [0, 0, 500, 500]);

Using the above methods will reduce your code length significantly. Additionally, you can also try the following ways to reduce the code length further: