Code covered by the BSD License  

Highlights from
MATLAB en zijn Symbolic Math Toolbox

image thumbnail

MATLAB en zijn Symbolic Math Toolbox

by

 

Leergang voor geïnteresseerden in de Computeralgebra - Systemen MATLAB en zijn Symbolic Math Toolbox

ContourIntegratie - kopie.m
%
% ***********************
% * ContourIntegratie.m *
% ***********************
%
% Deze m - file berekent oppervlakte - grootheden van
% polygoonvormige doorsneden. Dat wil zeggen :
%
% -    het oppervlak,
% -    de statisch momenten t.o.v. x - en y - as,
% -    de zwaartepuntsligging,
% -    de kwadratische oppervlakte - momenten t.o.v. het referentie 
%      assenstelsel,
% -    idem t.o.v. een parallel assenstelsel door het zwaarte punt.
%
% De zijden van het polygoon zijn genummerd : 1, 2, ..., N.
% De hoekpunten zijn genummerd :              1, 2, ..., N + 1.
%
% Zijde 1 heeft de hoekpunten 1 en 2, zijde 2 heeft 2 en 3, enz.. 
% In het algemeen geldt dat zijde i de hoekpunten i en i + 1 heeft.
%
% Omdat het een gesloten polygoon betreft, geldt :
%
%      x(N + 1) = x(1)
%      y(N + 1) = y(1)
%
% De invoer van dit script is vooraf geprepareerd in een file. Het  
% script vraagt om de naam hiervan, die vervolgens ingetoetst moet 
% worden, bijvoorbeeld : Contour.txt. (Dit is een string, zonder
% quotes.) In deze file staan, per hoekpunt, en in vrij formaat, de
% x - en y - coordinaten van het betreffende hoekpunt, dus 
% (bijvoorbeeld) als volgt :
%
%  0     0
% 12     0
% 12     4
%  4     4
%  4    16
% 16    16
% 16    20
%  0    20
%
% Vervolgens vraagt het script ook nog om een projectnaam. Ook dit
% is een string, die zonder quotes ingetoetst kan worden.
%

clear all; % workspace schoonmaken
clc;       % commandwindow schoonmaken
disp(' ');
filename = input(...
              'Naam van invoer - file, incl. extensie intoetsen : ', 's');
contname = input(...
              'Projectgegevens intoetsen                        : ', 's');
fid      = fopen(filename, 'r');      % openen van de opgegeven file
[a, Nco] = fscanf(fid, '%g %g', inf); % Nco succesvolle conversies, in de
                                      % een - dimensionale array a
N        = Nco / 2;                   % N is het aantal polygoonzijden 
for i = 1 : N
   x(i) = a(2 * i - 1);               % alle x - coordinaten uit a
   y(i) = a(2 * i);                   % alle y - coordinaten uit a
end % for i
x(N + 1) = x(1);
y(N + 1) = y(1);                      % hiermee is het polygoon gesloten

%
% Grafische uitvoer van het ingevoerde polygoon :
% -----------------------------------------------
%

maxx = max(x);
maxy = max(y);
minx = min(x);
miny = min(y); 
axis('square');
axis([minx, maxx, miny, maxy]);
axis('equal');
fill(x, y, 'g')
title(['Ingevoerd Contour : ', contname])
hold on

A   = 0;
Sx  = 0;
Sy  = 0;
Ix  = 0;
Iy  = 0;
Cxy = 0;

%
% Bepaling van de oppervlakte - grootheden in een loop :
% ------------------------------------------------------
%

for i = 1 : N
   d   = x(i) * y(i + 1) - x(i + 1) * y(i);
   A   = A   + d /  2;
   Sx  = Sx  + d /  6 * (y(i) + y(i + 1));
   Sy  = Sy  + d /  6 * (x(i) + x(i + 1));
   Ix  = Ix  + d / 12 * (y(i)^2 + y(i) * y(i + 1) + y(i + 1)^2);
   Iy  = Iy  + d / 12 * (x(i)^2 + x(i) * x(i + 1) + x(i + 1)^2);
   Cxy = Cxy + d / 24 * ...
             (x(i) * (2 * y(i) + y(i + 1)) + ...
              x(i + 1) * (y(i) + 2 * y(i + 1)));
end

%
% Bepaling van de zwaartepuntsligging :
% -------------------------------------
%

xZ  = Sy / A;
yZ  = Sx / A;
Ixx = Ix  - A * yZ^2;
Iyy = Iy  - A * xZ^2;
Ixy = Cxy - A * xZ * yZ;

%
% Uitvoer t.o.v. de referentie - assen en de zwaartepuntsligging :
% ----------------------------------------------------------------
%

disp(' ')
disp('Oppervlakte - grootheden t.o.v. de referentie - assen')
disp('-----------------------------------------------------')
disp(' ')
disp(sprintf('%s %15.2f', 'Oppervlak                    ', A   )) 
disp(sprintf('%s %15.2f', 'Statisch    Moment om x - as ', Sx  ))
disp(sprintf('%s %15.2f', 'Statisch    Moment om y - as ', Sy  ))
disp(sprintf('%s %15.2f', 'Kwadr. Opp. Moment om x - as ', Ix  ))
disp(sprintf('%s %15.2f', 'Kwadr. Opp. Moment om y - as ', Iy  ))
disp(sprintf('%s %15.2f', 'Centr. Opp. Moment om x - as ', Cxy ))
disp(' ')
disp('Ligging van het zwaartepunt')
disp('---------------------------')
disp(' ')
disp(sprintf('%s %15.2f', 'x - coordinaat zwaartepunt   ', xZ))
disp(sprintf('%s %15.2f', 'y - coordinaat zwaartepunt   ', yZ))

%
% Uitvoer t.o.v. assen door het zwaartepunt :
% -------------------------------------------
%

disp(' ')
disp('Oppervlakte - grootheden t.o.v. assen door het zwaartepunt')
disp('----------------------------------------------------------')
disp(' ')
disp(sprintf('%s %15.2f', 'Kwadr. Opp. Moment om x - as ', Ixx ))
disp(sprintf('%s %15.2f', 'Kwadr. Opp. Moment om y - as ', Iyy ))
disp(sprintf('%s %15.2f', 'Centr. Opp. Moment om x - as ', Ixy ))

%
% Berekening van de hoofdassen en de extreme waarden :
% ----------------------------------------------------
%
% De hoek alpha bepaalt de ligging van de maximale as
%

if Ixx ~= Iyy
   alpha = atan(-2 * Ixy / (Ixx - Iyy)) / 2;
   if Cxy ~= 0.0
      if sin(alpha * 2) / Ixy > 0
         alpha = alpha - sign(alpha) * pi / 2;
      end
   else 
      if Ixx < Iyy 
         alpha = pi / 2;
      end
   end 
else
   alpha = - sign(Ixy) * pi / 4;
end

%
% De hoek beta bepaalt de ligging van de minimale as
%

beta = alpha + pi / 2;

%
% Plotten van de hoofdassen
%

xMaxli = minx;
yMaxli = yZ - tan(alpha) * (xZ - minx);
if yMaxli < miny
   xMaxli = xZ - (yZ - miny) / tan(alpha);
   yMaxli = miny;
end;
xMaxre = maxx;
yMaxre = yZ + tan(alpha) * (maxx - xZ);
if yMaxre > maxy
   xMaxre = xZ + (maxy - yZ) / tan(alpha);
   yMaxre = maxy;
end;

xMinli = minx;
yMinli = yZ - tan(beta) * (xZ - minx);
if yMinli < miny
   xMinli = xZ - (yZ - miny) / tan(beta);
   yMinli = miny;
end;
xMinre = maxx;
yMinre = yZ + tan(beta) * (maxx - xZ);
if yMinre > maxy
   xMinre = xZ + (maxy - yZ) / tan(beta);
   yMinre = maxy;
end;

xMaxas = [xMaxli, xMaxre];
yMaxas = [yMaxli, yMaxre];
xMinas = [xMinli, xMinre];      
yMinas = [yMinli, yMinre];

plot(xMaxas, yMaxas, '--', xMinas, yMinas, ':')
legend('Profiel', 'maximale as ', 'minimale as ', -1)
plot(xZ, yZ, 'pentagram', 'Markersize', 12)
axis('equal');

Imax  = (Ixx + Iyy) / 2 + sqrt(((Ixx - Iyy) / 2)^2 + Ixy^2);
Imin  = (Ixx + Iyy) / 2 - sqrt(((Ixx - Iyy) / 2)^2 + Ixy^2);

disp(' ')
disp('Extreme waarden')
disp('---------------')
disp(' ')

disp(sprintf('%s %15.2f', 'Maximaal Kwadr. Opp. Moment  ', Imax))
disp(sprintf('%s %15.2f', 'Minimaal Kwadr. Opp. Moment  ', Imin))
disp(sprintf('%s %15.2f', 'Hoek alpha                   ', alpha))
disp(sprintf('%s %15.2f', 'Hoek beta                    ', beta ))

%
% *********************
% * ContourIntegratie *
% *********************
%




 

Contact us