Subs Function Memory Leak

Hi, I'm running a function (for discussion purposes I'll call it myfun.m) which utilizes both "solve()" and "subs()" function. Multiple runs of the myfun causes longer and longer execution times. I believe this issue is some sort of memory leak. If I close out of Matlab, open it back up, and execute myfun, it runs at the normal speed. If I continue executing it over and over it then begins to slow down. I was running 2007, but I'm upgrading to 2010. Has this problem been fixed? Is there some sort of clear function I can use that preserves my variables but clears whatever objects are being generated by subs or solve?

[Paulo Silva]

Please give us one simple example where the problem happens.

[Jan]

Please post any details. There is no known memory leak in SOLVE and SUBS. If there is a memory leak, it could be observed by inspecting the growing memory consumed Matlab. So you observe this?

Answer 1

[Walter Roberson]

solve() uses mupad which has its own memory area. Try using

reset(symengine)

[Walter Roberson]

Your code is somewhat long and complicated to go through.

If I were to go through it myself, I would be asking whether perhaps entire portions of it could be written symbolically and solve()'d, and then matlabFunction() used to convert that to MATLAB code that you could then just call with the actual numeric values.

Answer 2

[Michael Feeney]

In ver.2007 I was experiencing a memory leak. Now that I upgraded to ver.2010, there is no memory leak, however, the code takes longer to execute. I've attached the code and am certainly accepting of any suggestions to make it faster.

This is a pseudo ray tracer that projects a shadow onto a curved surface. The shadow is produced by a line with a circular or rectangular cross-section. The line is also defined by two points in space (p1 and p2).

    %%sblock.m
  % Description: This function is a spherical ray-tracer for circular or
  % rectangular cross-sectional elements on an arbitrary surface. The shadow 
  % area casted by arbitrarily oriented ciruclar cross-sectional elements 
  % within a surface's spherical reflection zone is calculated. Prime Focus 
  % must be located along [0 0 z]. Beam Elements that pass through the prime 
  % focus yeild an infinite blockage. WARNING: Do not have beam elements that
  % intersect the prime focus.
  %
  % |||||Date|||||||||||||||Author||||||||||||||||Changes|||||
  %   06-14-11       Michael E. Feeney         Original code.
  %
  % Input Variables:
  % p1 - Vertex of the element touching the surface of the primary [x1,y1,z1]
  % p2 - The other vertex of the element [x2,y2,z2]
  % d - Element Diameter

function [shadow_percentage] = sblock(p1,p2,w,d)
clear objects
clear functions
error(nargchk(3, 4, nargin))
if nargin == 3
    d = w;
end
if nargin == 4
    l=d;
end
t = 0:.25:1;   %Seed Size (Resolution)
len_t = length(t);          %Loop Length
syms u x_var y_var z_var;   %Parameterizing Variables

%%PRIMARY SURFACE AND RAY CONE DEFINITIONS
%Primary Surface Parameters
c1 = 0.05;               %Curvature
k = -1.00099;            %Konic Constant
D = 25;                  %Primary Surface Diameter
surface_area = 536.15;   %Primary Surface Area   
pf = [0 0 1/(2*c1)];     %Prime Focus
zmax = c1*((D/2)^2)/...  %Top Edge of the Primary Surface
    (1+sqrt(1-(1+k)*c1^2*((D/2)^2)));
surface = c1*(x_var^2+y_var^2)/(1+sqrt(1-(1+k)*c1^2*(x_var^2+y_var^2)))-z_var;
%Ray Cone Parameters
h1 = pf(3)-zmax;         %Height from Primary Diameter to Prime Focus
h2 = pf(3);              %Height from Primary Vertex to Prime Focus
r1 = D/2;                %Effective Radius 1
r2 = h2*r1/h1;           %Effective Radius 2
z_0 = pf(3);             %Ray Cone Parameter 1
c2 = r2/(pf(3));         %Ray Cone Parameter 2

%%CALCULATION OF BEAM INTERSECTION WITH PRIMARY SURFACE
%Checks to see the maximum height within the ray cone for p2(x,y)
height_check = c1*(p1(1)^2+p1(2)^2)/(1+sqrt(1-(1+k)*c1^2*(p1(1)^2+p1(2)^2)));
%Checks to see if the beam extends beyond the ray cone
if height_check > p1(3)
    slope1 = p2-p1;
    %Checks to see if the beam is vertical
    if (abs(slope1(1))<=1e-6 && abs(slope1(2))<=1e-6)
        vert = 1;
        x = p2(1);
        y = p2(2);
        z = p1(3) + u*slope1(3);
    else
        x = p1(1)+u*slope1(1);
        y = p1(2)+u*slope1(2);
        z = p1(3)+u*slope1(3);
    end
    %Calculates the intersection point
    prime_surface = subs(surface,x_var,x);
    prime_surface = subs(prime_surface,y_var,y);
    prime_surface = subs(prime_surface,z_var,z);
    uu = eval(solve(prime_surface,u));
    prime_point1 = [subs(x,u,uu(1)) subs(y,u,uu(1)) subs(z,u,uu(1))];
    prime_point2 = [subs(x,u,uu(2)) subs(y,u,uu(2)) subs(z,u,uu(2))];
    if (sqrt(prime_point1(1)^2+prime_point1(2)^2)-D/2<1e-6)
        p1 = prime_point1;
    else
        p1 = prime_point2;
    end

else
    slope1 = p2-p1;
    if (abs(slope1(1))<=1e-6 && abs(slope1(2))<=1e-6)
        vert = 1;
    end
    p1= p1;
end
fprintf('prime inter\n')
%%CALCULATION FOR A BEAM ALONG THE Z-AXIS
if (p1(1) == 0 && p1(2)==0 && p2(1)==0 && p2(2) ==0)
    if nargin == 4
        fprintf('DO NOT ALIGN RECTANGULAR BEAMS WITH PRIMARY SURFACE VERTEX AND PRIME FOCUS\n');
        shadow_percentage = sprintf('Invalid');
        return
    end
    %Checks to see if the beam goes through the prime focus
    if p2(3) >= pf(3)
        shadow_percentage = 100;
        return
    end
    %Determines line from prime focus to the top of the beam
    slope = (pf(3)-p2(3))/(pf(1)-(p2(1)+d/2));
    b = pf(3);
    z = slope*u+b;
    %Determines where the line intersects the primary surface
    prime_surface = subs(surface,x_var,x);
    prime_surface = subs(prime_surface,y_var,0);
    prime_surface = subs(prime_surface,z_var,z);
    %surface_inter =(0.05*(u^2)) / (1+sqrt(1-(1+k)*(u^2)*c1^2))-z;
    x = eval(solve(prime_surface,x_var));
    point1 = x(1);
    point2 = x(2);
    %Determines which point (if any) hits the bounded primary surface
    if (point1<=D/2)
        intersect_x = point1;
    elseif (point2<=D/2)
        intersect_x = point2;
    else
        shadow_percentage = 100;
        return
    end
    xs=intersect_x*(0:.1:10);
    %Calculates blockage percentage
    shadow_area =0;
    for i = 1:1:length(xs)-1
        shadow_area_fraction = pi*(xs(i+1)^2-xs(i)^2);
        shadow_area = shadow_area+shadow_area_fraction;
    end
    shadow_percentage = shadow_area/surface_area;
    return
end

%%CALCULATION FOR BEAM PORTION WITHIN THE RAY CONE
%Checks to see the maximum height within the ray cone for p2(x,y)
height_check = -sqrt((p2(1)^2+p2(2)^2)/c2^2)+z_0;
%Checks to see if the beam extends beyond the ray cone
if height_check < p2(3)
    slope1 = p2-p1;
    %Checks to see if the beam is vertical
    if (abs(slope1(1))<=1e-6 && abs(slope1(2))<=1e-6)
        vert = 1;
        x = p2(1);
        y = p2(2);
        z = p1(3) + u*slope1(3);
    else
        x = p1(1)+u*slope1(1);
        y = p1(2)+u*slope1(2);
        z = p1(3)+u*slope1(3);
    end
    %Calculates the intersection point
    cone_surface = (x^2+y^2)/c2^2-(z-z_0)^2;
    uu = eval(solve(cone_surface,u));
    cone_point1 = [subs(x,u,uu(1)) subs(y,u,uu(1)) subs(z,u,uu(1))];
    cone_point2 = [subs(x,u,uu(2)) subs(y,u,uu(2)) subs(z,u,uu(2))];
    if (sqrt(cone_point1(1)^2+cone_point1(2)^2)-D/2<1e-6)
        p2_act = cone_point1;
    else
        p2_act = cone_point2;
    end

else
    slope1 = p2-p1;
    if (abs(slope1(1))<=1e-6 && abs(slope1(2))<=1e-6)
        vert = 1;
    end
    p2_act = p2;
end
fprintf('ray cone\n')
%%CALCULATION OF THE SHADOW CENTER LINE
%Calculates appropriate slope for the beam
if exist('vert','var')==0
    beam_slope = p2_act-p1;
else
    beam_slope = [0 0 p2_act(3)-p1(3)];
end
%Calculates the coordinates for the beam within the ray cone
for i=1:1:len_t
    beam_center(i,:) = [p1(1)+beam_slope(1)*t(i), p1(2)+beam_slope(2)*t(i), p1(3)+beam_slope(3)*t(i)];
    %Vector Slope from Prime Focus to Element Point(m)
    slope = beam_center(i,:)-pf;
    %Parameterized Equations for Vector from Prime Focus to Beam
    x=pf(1)+u*slope(1);
    y=pf(2)+u*slope(2);
    z=pf(3)+u*slope(3);
    %Substitution of x,y,z into surface equation
    prime_surface = subs(surface,x_var,x);
    prime_surface = subs(prime_surface,y_var,y);
    prime_surface = subs(prime_surface,z_var,z);
    %Solves for Intersection Point Between Ray from Prime Focus to Surface.
    uu = eval(solve(prime_surface,u));
    if length(uu)==2
        surf_point1 = subs(z,u,uu(1));   %First point on surface
        surf_point2 = subs(z,u,uu(2));   %Second point on surface
        %Check to see which point is within the primary surface constraints
        if (norm(surf_point1,2)<norm(surf_point2,2))
            shadow_center(i,:) = [subs(x,u,uu(1)) subs(y,u,uu(1)) surf_point1];
        else
            shadow_center(i,:) = [subs(x,u,uu(2)) subs(y,u,uu(2)) surf_point2];
        end
    else
        surf_point1 = subs(z,u,uu(1));   %First point on surface
        shadow_center(i,:) = [subs(x,u,uu(1)) subs(y,u,uu(1)) surf_point1];
    end
end
fprintf('center line\n')
%%TRANSFORMATION MATRIX (LOCAL TO GLOBAL)
%Length of Element
L = norm(p2_act-p1,2);
%Local-Global Coordinate Transformation
if exist('vert','var')
    T = [0 0 1; 0 1 0; -1 0 0; ];
else
    l_x = (p2_act(1)-p1(1))/L;
    m_x = (p2_act(2)-p1(2))/L;
    n_x = (p2_act(3)-p1(3))/L;
    D = (l_x^2+m_x^2)^0.5;
    l_y = -m_x/D;
    m_y = l_x/D;
    n_y = 0;
    l_z = -l_x*n_x/D;
    m_z = -m_x*n_x/D;
    n_z = D;
    T = [l_x m_x n_x; l_y m_y n_y; l_z m_z n_z];
end

%%CALCULATION OF SIDE 1 EDGE SHADOW (CIRCULAR CROSS-SECTION)
if nargin == 3
    %Phi Iteration Mid-Beam (Side1)
    r = d/2;
    phi = 0:pi/30:(2*pi-pi/30);
    slope_check =0;
    len_2 = ceil(length(t)/2);
    fprintf('phi iter begin\n')
    for j = len_2
        for i = 1:1:length(phi)
            x = t(j)*L;
            y = r*sin(phi(i));
            z = r*cos(phi(i));
            p3 = [x y z];
            check = T'*p3'+p1';
            element_edge= T'*p3'+p1';
            slope3 = element_edge'-pf;
            x1=pf(1)+u*slope3(1);
            y1=pf(2)+u*slope3(2);
            z1=pf(3)+u*slope3(3);
            shadow = c1*(x1^2+y1^2)/(1+sqrt(1-(1+k)*c1^2*(x1^2+y1^2)))-z1;
            uu1 = eval(solve(shadow,u));
            if length(uu1) == 2
                shadow1 = [subs(x1,u,uu1(1)) subs(y1,u,uu1(1)) subs(z1,u,uu1(1))];
                shadow2 = [subs(x1,u,uu1(2)) subs(y1,u,uu1(2)) subs(z1,u,uu1(2))];
                if (shadow1(3)<shadow2(3))
                    shadow_point1 = shadow1;
                elseif (shadow2(3)<shadow1(3))
                    shadow_point1 = shadow2;
                end
            else
                shadow1 = [subs(x1,u,uu1(1)) subs(y1,u,uu1(1)) subs(z1,u,uu1(1))];
                shadow_point1 = shadow1;
            end
            if j == len_2
                slope_check1 = (shadow_center(1,2)-shadow_point1(2))/(shadow_center(1,1)-shadow_point1(1));
                slope_check2 = (shadow_center(1,2)-shadow_center(2,2))/(shadow_center(1,1)-shadow_center(2,1));
                slope_diff = abs(slope_check1-slope_check2);
                store(i,:) = [slope_diff phi(i)];
                if slope_diff>slope_check
                    slope_check = slope_diff;
                    phi2 = phi(i);
                    shadow_edge1(j,:) = shadow_point1;
                end
            end
        end
    end
    fprintf('phi iter done\n')
    %Calculation for the rest of Side 1
    for i=1:1:len_t
        if i~=len_2
            x = t(i)*L;
            y = r*sin(phi2);
            z = r*cos(phi2);
            p3 = [x y z];
            element_edge= T'*p3'+p1';
            slope3 = element_edge'-pf;
            x1=pf(1)+u*slope3(1);
            y1=pf(2)+u*slope3(2);
            z1=pf(3)+u*slope3(3);
            shadow = c1*(x1^2+y1^2)/(1+sqrt(1-(1+k)*c1^2*(x1^2+y1^2)))-z1;
            uu1 = eval(solve(shadow,u));
            if length(uu1) == 2
                shadow1 = [subs(x1,u,uu1(1)) subs(y1,u,uu1(1)) subs(z1,u,uu1(1))];
                shadow2 = [subs(x1,u,uu1(2)) subs(y1,u,uu1(2)) subs(z1,u,uu1(2))];
                if (shadow1(3)<shadow2(3))
                    shadow_point1 = shadow1;
                elseif (shadow2(3)<shadow1(3))
                    shadow_point1 = shadow2;
                end
            else
                shadow1 = [subs(x1,u,uu1(1)) subs(y1,u,uu1(1)) subs(z1,u,uu1(1))];
                shadow_point1 = shadow1;
            end
            shadow_edge1(i,:) = shadow_point1;
        end
    end
    [~, locs] = findpeaks(store(:,1));
    if phi(locs(1)) == phi2
        phi2 = phi(locs(2));
    else
        phi2 = phi(locs(1));
    end
    fprintf('side done\n')
    %Calculation for Side 2
    for i=1:1:len_t
        x = t(i)*L;
        y = r*sin(phi2);
        z = r*cos(phi2);
        p3 = [x y z];
        element_edge= T'*p3'+p1';
        slope3 = element_edge'-pf;
        x1=pf(1)+u*slope3(1);
        y1=pf(2)+u*slope3(2);
        z1=pf(3)+u*slope3(3);
        shadow = c1*(x1^2+y1^2)/(1+sqrt(1-(1+k)*c1^2*(x1^2+y1^2)))-z1;
        uu1 = eval(solve(shadow,u));
        if length(uu1) == 2
            shadow1 = [subs(x1,u,uu1(1)) subs(y1,u,uu1(1)) subs(z1,u,uu1(1))];
            shadow2 = [subs(x1,u,uu1(2)) subs(y1,u,uu1(2)) subs(z1,u,uu1(2))];
            if (shadow1(3)<shadow2(3))
                shadow_point1 = shadow1;
            elseif (shadow2(3)<shadow1(3))
                shadow_point1 = shadow2;
            end
        else
            shadow1 = [subs(x1,u,uu1(1)) subs(y1,u,uu1(1)) subs(z1,u,uu1(1))];
            shadow_point1 = shadow1;
        end
        shadow_edge2(i,:) = shadow_point1;
    end
end
fprintf('side done\n')
%%CALCULATION OF SIDE 1 EDGE SHADOW (RECTANGULAR CROSS-SECTION)
if nargin == 4
    len_2 = ceil(length(t)/2);
    x = t(len_2)*L;
    w = w/2;
    l = l/2;
    loc_corner1 = [x w l];
    loc_corner2 = [x -w l];
    glo_corner1 = T'*loc_corner1'+p1';
    glo_corner2 = T'*loc_corner2'+p1';
    glo_corners = [glo_corner1'; glo_corner2'];
    for i = 1:1:2
        slope3 = glo_corners(i,:)-pf;
        x1=pf(1)+u*slope3(1);
        y1=pf(2)+u*slope3(2);
        z1=pf(3)+u*slope3(3);
        shadow = c1*(x1^2+y1^2)/(1+sqrt(1-(1+k)*c1^2*(x1^2+y1^2)))-z1;
        uu1 = eval(solve(shadow,u));
        if length(uu1) == 2
            shadow1 = [subs(x1,u,uu1(1)) subs(y1,u,uu1(1)) subs(z1,u,uu1(1))];
            shadow2 = [subs(x1,u,uu1(2)) subs(y1,u,uu1(2)) subs(z1,u,uu1(2))];
            if (shadow1(3)<shadow2(3))
                shadow_point1 = shadow1;
            elseif (shadow2(3)<shadow1(3))
                shadow_point1 = shadow2;
            end
        else
            shadow1 = [subs(x1,u,uu1(1)) subs(y1,u,uu1(1)) subs(z1,u,uu1(1))];
            shadow_point1 = shadow1;
        end
        shadow_edge(i,:) = shadow_point1;
    end
    slope_check1 = (shadow_center(1,2)-shadow_edge(1,2))/(shadow_center(1,1)-shadow_edge(1,1));
    slope_check2 = (shadow_center(1,2)-shadow_edge(2,2))/(shadow_center(1,1)-shadow_edge(2,1));
    slope_check3 = (shadow_center(1,2)-shadow_center(len_2,2))/(shadow_center(1,1)-shadow_center(len_2,1));
    slope_diff1 = abs(slope_check1-slope_check3);
    slope_diff2 = abs(slope_check2-slope_check3);
    check_inline1 = abs(abs(glo_corner1(1))-abs(glo_corner2(1)));
    check_inline2 = abs(abs(glo_corner1(2))-abs(glo_corner2(2)));
    check_inline = check_inline1+check_inline2;
    if (slope_diff1>slope_diff2 || (check_inline<1e-10))
        shadow_edge1(2,:) = shadow_edge(1,:);
        for i=1:1:len_t
            x = t(i)*L;
            loc_corner1 = [x w l];
            loc_corner2 = [x -w l];
            loc_corner3 = [x -w -l];
            loc_corner4 = [x w -l];
            glo_corner1 = T'*loc_corner1'+p1';
            glo_corner2 = T'*loc_corner2'+p1';
            glo_corner3 = T'*loc_corner3'+p1';
            glo_corner4 = T'*loc_corner4'+p1';
            if (((check_inline)<1e-10) && i~=len_t)
                slope1 = glo_corner4'-pf;
                slope3 = glo_corner3'-pf;
            elseif ((check_inline<1e-10) && i==len_t)
                slope1 = glo_corner1'-pf;
                slope3 = glo_corner2'-pf;
            elseif (check_inline>1e-10)
                slope1 = glo_corner1'-pf;
                slope3 = glo_corner3'-pf;
            end
            x1=pf(1)+u*slope1(1);
            y1=pf(2)+u*slope1(2);
            z1=pf(3)+u*slope1(3);
            x3=pf(1)+u*slope3(1);
            y3=pf(2)+u*slope3(2);
            z3=pf(3)+u*slope3(3);
            shadow_1 = c1*(x1^2+y1^2)/(1+sqrt(1-(1+k)*c1^2*(x1^2+y1^2)))-z1;
            shadow_3 = c1*(x3^2+y3^2)/(1+sqrt(1-(1+k)*c1^2*(x3^2+y3^2)))-z3;
            uu1 = eval(solve(shadow_1,u));
            uu3 = eval(solve(shadow_3,u));
            if length(uu1) == 2
                shadow1 = [subs(x1,u,uu1(1)) subs(y1,u,uu1(1)) subs(z1,u,uu1(1))];
                shadow2 = [subs(x1,u,uu1(2)) subs(y1,u,uu1(2)) subs(z1,u,uu1(2))];
                if (shadow1(3)<shadow2(3))
                    shadow_point1 = shadow1;
                elseif (shadow2(3)<shadow1(3))
                    shadow_point1 = shadow2;
                end
            else
                shadow1 = [subs(x1,u,uu1(1)) subs(y1,u,uu1(1)) subs(z1,u,uu1(1))];
                shadow_point1 = shadow1;
            end
            if length(uu3) == 2
                shadow1 = [subs(x3,u,uu3(1)) subs(y3,u,uu3(1)) subs(z3,u,uu3(1))];
                shadow2 = [subs(x3,u,uu3(2)) subs(y3,u,uu3(2)) subs(z3,u,uu3(2))];
                if (shadow1(3)<shadow2(3))
                    shadow_point3 = shadow1;
                elseif (shadow2(3)<shadow1(3))
                    shadow_point3 = shadow2;
                end
            else
                shadow1 = [subs(x3,u,uu3(1)) subs(y3,u,uu3(1)) subs(z3,u,uu3(1))];
                shadow_point3 = shadow1;
            end
            shadow_edge1(i,:) = shadow_point1;
            shadow_edge2(i,:) = shadow_point3;
        end
    else
        for i=1:1:len_t
            x = t(i)*L;
            loc_corner2 = [x -w l];
            glo_corner2 = T'*loc_corner2'+p1';
            loc_corner4 = [x w -l];
            glo_corner4 = T'*loc_corner4'+p1';
            slope2 = glo_corner2'-pf;
                slope4 = glo_corner4'-pf;
                x2=pf(1)+u*slope2(1);
                y2=pf(2)+u*slope2(2);
                z2=pf(3)+u*slope2(3);
                x4=pf(1)+u*slope4(1);
                y4=pf(2)+u*slope4(2);
                z4=pf(3)+u*slope4(3);
                shadow_2 = c1*(x2^2+y2^2)/(1+sqrt(1-(1+k)*c1^2*(x2^2+y2^2)))-z2;
                shadow_4 = c1*(x4^2+y4^2)/(1+sqrt(1-(1+k)*c1^2*(x4^2+y4^2)))-z4;
                uu2 = eval(solve(shadow_2,u));
                uu4 = eval(solve(shadow_4,u));
                if length(uu2) == 2
                    shadow1 = [subs(x2,u,uu2(1)) subs(y2,u,uu2(1)) subs(z2,u,uu2(1))];
                    shadow2 = [subs(x2,u,uu2(2)) subs(y2,u,uu2(2)) subs(z2,u,uu2(2))];
                    if (shadow1(3)<shadow2(3))
                        shadow_point2 = shadow1;
                    elseif (shadow2(3)<shadow1(3))
                        shadow_point2 = shadow2;
                    end
                else
                    shadow1 = [subs(x2,u,uu2(1)) subs(y2,u,uu2(1)) subs(z2,u,uu2(1))];
                    shadow_point2 = shadow1;
                end
                if length(uu4) == 2
                    shadow1 = [subs(x4,u,uu4(1)) subs(y4,u,uu4(1)) subs(z4,u,uu4(1))];
                    shadow2 = [subs(x4,u,uu4(2)) subs(y4,u,uu4(2)) subs(z4,u,uu4(2))];
                    if (shadow1(3)<shadow2(3))
                        shadow_point4 = shadow1;
                    elseif (shadow2(3)<shadow1(3))
                        shadow_point4 = shadow2;
                    end
                else
                    shadow1 = [subs(x4,u,uu4(1)) subs(y4,u,uu4(1)) subs(z4,u,uu4(1))];
                    shadow_point4 = shadow1;
                end
                shadow_edge1(i,:) = shadow_point2;
                shadow_edge2(i,:) = shadow_point4;
        end
    end
end

%%AREA CALCULATION
shadow_area = 0;
for i=1:1:len_t-1
    avg_width1 = sqrt(dot(shadow_edge1(i,1:3)-shadow_edge2(i,:),shadow_edge1(i,1:3)-shadow_edge2(i,:)));
    avg_width2 = sqrt(dot(shadow_edge1(i+1,1:3)-shadow_edge2(i+1,:),shadow_edge1(i+1,1:3)-shadow_edge2(i+1,:)));
    avg_length1 = sqrt(dot(shadow_edge1(i,1:3)-shadow_edge1(i+1,1:3),shadow_edge1(i,1:3)-shadow_edge1(i+1,1:3)));
    avg_length2 = sqrt(dot(shadow_edge2(i,:)-shadow_edge2(i+1,:),shadow_edge2(i,:)-shadow_edge2(i+1,:)));
    avg_width = (avg_width1+avg_width2)/2;
    avg_length = (avg_length1+avg_length2)/2;

      area(i) = avg_width*avg_length;
      shadow_area = shadow_area+area(i);
  end
  shadow_percentage = 100*(shadow_area / surface_area);
  end

[Walter Roberson]

Style note:

eval(solve(...))

should be replaced with

double(subs(solve(...)))

or

double(solve(...)))

depending on whether the solve() returns symbolic numbers only or formula with names whose current value you want substituted in.

[Michael Feeney]

Hey Walter,

Thanks for your input. I'm not sure I fully understand why answer_var = eval(solve(eqn,u)) is wrong. What this combination does is the following:

solve - solves a second order equation 'eqn' for the two values of 'u'

eval - puts these two values of 'u' into a vector 'answer_var'

in other words, I need the two values of 'u' to be stored into a vector, rather than echoed out.

[Walter Roberson]

eval() invokes the full matlab parser and JIT and threaded interpreter, and is thus a heavy-weight method of dealing with the situation.

subs(solve(...)) will look for variables in the result of solve() and will, if possible, replace them with values from the current workspace, and then will do whatever further calculations are practical. That should leave you with a symbolic number or a symbolic formula involving only constants; double() will then convert that symbolic number or constant symbolic formula in to a MATLAB-level double-precision number. The symbolic engine will handle this efficiently.

[Michael Feeney]

Wow. Okay, that actually makes quite a bit of sense. I was unaware of the computational intensity of eval. Looks like I'll do a Crtl+F and replace eval(solve( with double(subs(solve(. Once again, I appreciate your input, I hope this reduces the execution time. I'll let you know.

Thanks.