Warning: Error updating FunctionLine in using fplot

Hi all, I wrote a function, using the PDE modeler app, that takes a radius as input and solves a PDE on a shape depending on r, and integrates the solution on the same shape. The function by itself seems to be working. However, when I try tu use fplot to plot it, it gives me the following warning: "Warning: Error updating FunctionLine. The following error was reported evaluating the function in FunctionLine update: Dimensions of arrays being concatenated are not consistent.", and I'm not really sure why that is. The code is below.

Thanks in advance!

fplot(@(r) Ttr(r), [1-1/sqrt(3), 0.47])

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

function [rt] = Ttr(r)
    C1=[1; 1-r; 0; r; 0; 0; 0; 0];
    P1=[2; 3; 0; (1-2*r)/(1-r); (1-2*r)/(1-r); 0; r/(1-r)*sqrt(1-2*r); -r/(1-r)*sqrt(1-2*r)];
    dg=[C1, P1];
    ns = char('C1','P1');
    ns = ns';
    sf='C1+P1';
    [dl,bt] = decsg(dg,sf,ns);
    [dl1,~] = csgdel(dl,bt);
    %pdegplot(dl1,"EdgeLabels","on","FaceLabels","on")
    %e=input("Inserire il numero di lati della figura: ");
     
    model = createpde();
    geometryFromEdges(model,dl1);
    applyBoundaryCondition(model,"dirichlet","Edge", 1:6,"u",0);
    specifyCoefficients(model,"m",0,"d",0,"c",1,"a",0,"f",1);
    mesh=generateMesh(model, Hmax=0.05, Hmin=0.00005, GeometricOrder='linear');
    [p,~,t] = meshToPet(mesh);
    u=solvepde(model);
    
    rt=0;
    k=length(t(1,:));
    for j=1:k
        A=p(:,t(1,j));
        B=p(:,t(2,j));
        C=p(:,t(3,j));
        M=[A(1),A(2),1; B(1),B(2),1; C(1),C(2),1];
        rt=rt+1/2*(u.NodalSolution(t(1,j))+u.NodalSolution(t(2,j))+u.NodalSolution(t(3,j)))/3 * abs(det(M));
    end
end

Answer 1

[dpb]

function [rt] = Ttr(r)
    C1=[1; 1-r; 0; r; 0; 0; 0; 0];
    P1=[2; 3; 0; (1-2*r)/(1-r); (1-2*r)/(1-r); 0; r/(1-r)*sqrt(1-2*r); -r/(1-r)*sqrt(1-2*r)];
    dg=[C1, P1];
    ns = char('C1','P1');
    ns = ns';
    sf='C1+P1';
    [dl,bt] = decsg(dg,sf,ns);
    [dl1,~] = csgdel(dl,bt);
    %pdegplot(dl1,"EdgeLabels","on","FaceLabels","on")
    %e=input("Inserire il numero di lati della figura: ");
     
    model = createpde();
    geometryFromEdges(model,dl1);
    applyBoundaryCondition(model,"dirichlet","Edge", 1:6,"u",0);
    specifyCoefficients(model,"m",0,"d",0,"c",1,"a",0,"f",1);
    mesh=generateMesh(model, Hmax=0.05, Hmin=0.00005, GeometricOrder='linear');
    [p,~,t] = meshToPet(mesh);
    u=solvepde(model);
    
    rt=0;
    k=length(t(1,:));
    for j=1:k
        A=p(:,t(1,j));
        B=p(:,t(2,j));
        C=p(:,t(3,j));
        M=[A(1),A(2),1; B(1),B(2),1; C(1),C(2),1];
        rt=rt+1/2*(u.NodalSolution(t(1,j))+u.NodalSolution(t(2,j))+u.NodalSolution(t(3,j)))/3 * abs(det(M));
    end
end
%fplot(@(r) Ttr(r), [1-1/sqrt(3), 0.47])

Let's ilustrate the problem...

r=[1-1/sqrt(3), 0.47];
rt = Ttr(r);

Error using vertcat
Dimensions of arrays being concatenated are not consistent.

Error in solution>Ttr (line 2)
C1=[1; 1-r; 0; r; 0; 0; 0; 0];

Now you can see it clearly; the issue is that fplot tries to evaluate the functional as being vectorized and so passes more than one argument at a time to the function. But, your function is constructed such that it expects (and requires) only a single value of r each time it is called; it tries/needs to build a vector of constants based on that input value, but if that value is a vector instead of a single value, it will be unable to do so.

You could bypass that particular error if you passed a column vector instead of a row vector, but the function still wouldn't work because then the definition of C1, P1 wouldn't match the form intended, even if it would run (which I didn't try, but I suspect wouldn't).

What you would need to do would be to test inside Ttr() for a vector input and then loop over those values individually and return the result of each in an array of the same dimensions as r going in.

The first edit forgot to then only pick the one r value from the vector of values passed each loop so while it checked and allocated, it still tried to catenate the whole input array r.

The following corrected code renames the r input argument to rin and allocates the output array based on its size then sets r to the single value for each pass through the loop.

function [rt] = Ttr(rin)                        % change input argument name so can keep r later
  assert(isvector(rin),'Inputs must be vector') % abort, tell user if not vector input
  N=numel(rin);                                 % find out how many are requested
  rt=nan(size(rin));                            % preallocate; nan will let know if failed
  for i=1:N
    r=rin(i);                                   % pick up the one element for each pass
    K1=(1-2*r)/(1-r);                           % compute repeated constants
    K2=r/(1-r)*sqrt(1-2*r);                     % may be minute efficiency
    C1=[1; 1-r; 0; r; 0; 0; 0; 0];
    P1=[2; 3; 0; K1; K1; 0; K2; -K2];
    dg=[C1, P1];
    ns = char('C1','P1');
    ns = ns';
    sf='C1+P1';
    [dl,bt] = decsg(dg,sf,ns);
    [dl1,~] = csgdel(dl,bt);
    %pdegplot(dl1,"EdgeLabels","on","FaceLabels","on")
    %e=input("Inserire il numero di lati della figura: ");
     
    model = createpde();
    geometryFromEdges(model,dl1);
    applyBoundaryCondition(model,"dirichlet","Edge", 1:6,"u",0);
    specifyCoefficients(model,"m",0,"d",0,"c",1,"a",0,"f",1);
    mesh=generateMesh(model, Hmax=0.05, Hmin=0.00005, GeometricOrder='linear');
    [p,~,t] = meshToPet(mesh);
    u=solvepde(model);
    
    rt(i)=0;                % initialize the accumulator
    k=length(t(1,:));
    for j=1:k
        A=p(:,t(1,j));
        B=p(:,t(2,j));
        C=p(:,t(3,j));
        M=[A(1),A(2),1; B(1),B(2),1; C(1),C(2),1];
        rt(i)=rt(i)+1/2*(u.NodalSolution(t(1,j))+u.NodalSolution(t(2,j))+u.NodalSolution(t(3,j)))/3 * abs(det(M));
    end
  end  % loop for number of r input...    
end

[Simone]

Thank you for your answer. I think I understand your point, but I have some doubts: the first is the functioning of fplot. I read its help page, and with the syntax fplot(f,xinterval) I assumed it would have been a pointwise calculation in a discretization of the interval, but I'm trusting you on how it works (also given the fact that my code does not work :) ).

Then, in your code in line 4, is rt=nan(size(rt)); correct or should it have been rt=nan(size(r));?

Either way, the code still doesn't seem to be working (unless I'm making some other mistake), giving the same error.

[dpb]

Yes, my bad on using rt instead of r for the preallocation size...fixed above.

While it's not at the top level in the Description or Syntax sections (where it probably should be), the description of the function input argument says this...

..."The function must accept a vector input argument and return a vector output argument of the same size. Use array operators instead of matrix operators for the best performance. For example, use .* (times) instead of * (mtimes)."

The reason you got the same error is that I also forgot to only address the individual r elements so it was still trying to do the same thing in the creation of C,P...to make less editing, I'll rename the input argument so can still use r internally...

Let's try it again now with the above correction...

function [rt] = Ttr(rin)
  assert(isvector(rin),'Inputs must be vector')     % abort, tell user if not vector input
  N=numel(rin);                                     % find out how many
  rt=nan(size(rin));                                % preallocate; nan will let know if failed
  for i=1:N
    r=rin(i);                                       % get the one element wanted each pass
    C1=[1; 1-r; 0; r; 0; 0; 0; 0];
    P1=[2; 3; 0; (1-2*r)/(1-r); (1-2*r)/(1-r); 0; r/(1-r)*sqrt(1-2*r); -r/(1-r)*sqrt(1-2*r)];
    dg=[C1, P1];
    ns = char('C1','P1');
    ns = ns';
    sf='C1+P1';
    [dl,bt] = decsg(dg,sf,ns);
    [dl1,~] = csgdel(dl,bt);
    %pdegplot(dl1,"EdgeLabels","on","FaceLabels","on")
    %e=input("Inserire il numero di lati della figura: ");
     
    model = createpde();
    geometryFromEdges(model,dl1);
    applyBoundaryCondition(model,"dirichlet","Edge", 1:6,"u",0);
    specifyCoefficients(model,"m",0,"d",0,"c",1,"a",0,"f",1);
    mesh=generateMesh(model, Hmax=0.05, Hmin=0.00005, GeometricOrder='linear');
    [p,~,t] = meshToPet(mesh);
    u=solvepde(model);
    
    rt(i)=0;                % initialize the accumulator
    k=length(t(1,:));
    for j=1:k
        A=p(:,t(1,j));
        B=p(:,t(2,j));
        C=p(:,t(3,j));
        M=[A(1),A(2),1; B(1),B(2),1; C(1),C(2),1];
        rt(i)=rt(i)+1/2*(u.NodalSolution(t(1,j))+u.NodalSolution(t(2,j))+u.NodalSolution(t(3,j)))/3 * abs(det(M));
    end
  end  % loop for number of r input...    
end
r=[1-1/sqrt(3), 0.47];
Ttr(r)

ans = 1×2

0.0135 0.0195

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

and voila! now no errors.

The alternative solution would be to leave the original function as is and the write a wrapper function around it that does the vectorization that is the function called in fplot.

I'll make the correction above, but leave this here for the time being for clarification. I knew we needed to do that, just got ahead of myself...

[Simone]

Now it works, thank you!

[Torsten]

The easiest way is to use

fplot(@(r) arrayfun(@(r)Ttr(r),r), [1-1/sqrt(3), 0.47])

This will supply the r-array from "fplot" elementwise to your function "Ttr".

Answer 2

[Torsten]

In R2024b, I get the warning (see above):

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

but fplot nevertheless works. The reason for this warning is that "fplot" calls your function with an array of values for r as input, but your function is not able to handle this array input.