Error when multiplying an evaluated function

I am trying to create a matrix K from smaller matricies Kelem. When I make a function call to Kelem I get the following error. I know it has something to do with the function A that represents different values depending on the parameter x. Can someone help me resolve this?

Error using integralCalc/finalInputChecks

Input function must return 'double' or 'single' values. Found 'function_handle'.

Error in integralCalc/iterateScalarValued (line 315)

finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 132)

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)

[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)

Q = integralCalc(fun,a,b,opstruct);

Error in hw3_872>Kelem (line 170)

out = 1/2/Lelem*[integral(@(eta) A(x),-1,1) integral(@(eta) -A(x),-1,1); integral(@(eta) -A(x),-1,1) integral(@(eta) A(x),-1,1)];

Error in hw3_872 (line 17)

K(elem(e,1):elem(e,2), elem(e,1):elem(e,2)) = K(elem(e,1):elem(e,2), elem(e,1):elem(e,2)) + Kelem(e/N,N,L);

N = 100;
L = 1;
node = linspace(0,L,N+1);
nodes = linspace(0,L,1000*(N+1));
elem = zeros(N,2);
for e = 1:N
    elem(e,1) = e;
    elem(e,2) = e+1;
end

K = sparse(N+1,N+1);
R = zeros(N+1,1);
k = 32;
L = 1;

for e = 1:N
    K(elem(e,1):elem(e,2), elem(e,1):elem(e,2)) = K(elem(e,1):elem(e,2), elem(e,1):elem(e,2)) + Kelem(e/N,N,L);
    R(elem(e,1):elem(e,2)) = R(elem(e,1):elem(e,2)) + Relem((elem(e,1)-1).*(L/N), (elem(e,2)-1).*(L/N),N,k,L);
end
function out = Kelem(x,N,L)
    Lelem = L/N;
    out = 1/2/Lelem*[integral(@(eta) A(x),-1,1) integral(@(eta) -A(x),-1,1); integral(@(eta) -A(x),-1,1) integral(@(eta) A(x),-1,1)];
end

function out = Relem(x1,x2,N,k,L)
    Lelem = L/N;
    phi1 = @(eta) 1/2*(1-eta);
    phi2 = @(eta) 1/2*(1+eta);
    out = [integral(@(eta) phi1(eta).*fun(eta,x1,x2,k,L), -1, 1).*(Lelem/2.0); integral(@(eta) phi2(eta).*fun(eta,x1,x2,k,L), -1, 1).*(Lelem/2.0)];
end

function out = fun(eta,x1,x2,k,L)
    out = -k^2*cos(2*pi*k*parametric(eta,x1,x2)/L);
end

function out = parametric(eta,x1,x2)
    out = (1-eta)/2.0*x1 + (1+eta)/2.0*x2;
end

function out = A(x)
    if x < .1
        A = @(x) 2 + 0*fun(eta,x1,x2,k,L);
    elseif x < .2
        A = @(x) 1.5 + 0*fun(eta,x1,x2,k,L);
    elseif x < .3
        A = @(x) 1.75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .4
        A = @(x) 1.5 + 0*fun(eta,x1,x2,k,L);
    elseif x < .5
        A = @(x) 3.75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .6
        A = @(x) .75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .7
        A = @(x) 1.25 + 0*fun(eta,x1,x2,k,L);
    elseif x < .8
        A = @(x) .75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .9
        A = @(x) 2 + 0*fun(eta,x1,x2,k,L);
    elseif x <= 1
        A = @(x) 1 + 0*fun(eta,x1,x2,k,L);
    end
    out = A;
end

[Stephen23]

The basic problem is that you have not parameterized your functions correctly:

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

As your code is written now, it makes no sense to call INTEGRAL with a function handle that competely ignores its input argument. You need to use its input argument and parameterize all of your functions handles correctly.

[Tyler Michael Paparella]

Thank you for the response. I read the article and made some changes to my code, but I am still getting the same error. I am not sure what you mean by A(x) ignoring its input argument. x determines the constant function A. Maybe you mean I am multiplying fun() by 0, but even if I multiply by 1, for example, I get the same error.

N = 100;
L = 1;
node = linspace(0,L,N+1);
nodes = linspace(0,L,1000*(N+1));
elem = zeros(N,2);
for e = 1:N
    elem(e,1) = e;
    elem(e,2) = e+1;
end

K = sparse(N+1,N+1);
R = zeros(N+1,1);
k = 32;
L = 1;

for e = 1:N
    K(elem(e,1):elem(e,2), elem(e,1):elem(e,2)) = K(elem(e,1):elem(e,2), elem(e,1):elem(e,2)) + Kelem(e/N,(elem(e,1)-1).*(L/N), (elem(e,2)-1).*(L/N),k,N,L);
    R(elem(e,1):elem(e,2)) = R(elem(e,1):elem(e,2)) + Relem((elem(e,1)-1).*(L/N), (elem(e,2)-1).*(L/N),N,k,L);
end

function out = Kelem(x,x1,x2,k,N,L)
    Lelem = L/N;
    out = 1/2/Lelem*[integral(@(eta) A(x,x1,x2,k,L),-1,1) integral(@(eta) -A(x,x1,x2,k,L),-1,1); integral(@(eta) -A(x,x1,x2,k,L),-1,1) integral(@(eta) A(x,x1,x2,k,L),-1,1)];
end

function out = Relem(x1,x2,N,k,L)
    Lelem = L/N;
    phi1 = @(eta) 1/2*(1-eta);
    phi2 = @(eta) 1/2*(1+eta);
    out = [integral(@(eta) phi1(eta).*fun(eta,x1,x2,k,L), -1, 1).*(Lelem/2.0); integral(@(eta) phi2(eta).*fun(eta,x1,x2,k,L), -1, 1).*(Lelem/2.0)];
end

function out = fun(eta,x1,x2,k,L)
    out = -k^2*cos(2*pi*k*parametric(eta,x1,x2)/L);
end

function out = parametric(eta,x1,x2)
    out = (1-eta)/2.0*x1 + (1+eta)/2.0*x2;
end

function out = A(x,x1,x2,k,L)
    if x < .1
        out = @(eta) 2 + 0*fun(eta,x1,x2,k,L);
    elseif x < .2
        out = @(eta) 1.5 + 0*fun(eta,x1,x2,k,L);
    elseif x < .3
        out = @(eta) 1.75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .4
        out = @(eta) 1.5 + 0*fun(eta,x1,x2,k,L);
    elseif x < .5
        out = @(eta) 3.75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .6
        out = @(eta) .75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .7
        out = @(eta) 1.25 + 0*fun(eta,x1,x2,k,L);
    elseif x < .8
        out = @(eta) .75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .9
        out = @(eta) 2 + 0*fun(eta,x1,x2,k,L);
    elseif x <= 1
        out = @(eta) 1 + 0*fun(eta,x1,x2,k,L);
    end
end

[Stephen23]

" I am not sure what you mean by A(x) ignoring its input argument. x determines the constant function A. Maybe you mean I am multiplying fun() by 0, but even if I multiply by 1, for example, I get the same error."

I did not write anything about "A(x) ignoring its input argument", nor is it clear to me what you mean by that.

You define lots of anonymous functions, e.g.:

@(eta) A(x,x1,x2,k,L)

Exactly as I wrote: the input ETA is totally unused, so there is no point in calling INTEGRAL with those function handles.

Lets take a look at an example from the INTEGRAL documentation. This is the very first example:

@(x) exp(-x.^2).*log(x).^2;

Note how X is an input argument and how it is used within the anonymous function. Compare against your code.

While there are certainly circumstances where it can be useful to define an anonymous function with input arguments that get ignored, INTEGRAL is not one of them.

Answer 1

[VBBV]

N = 100;
L = 1;
node = linspace(0,L,N+1);
nodes = linspace(0,L,1000*(N+1));
elem = zeros(N,2);
for e = 1:N
    elem(e,1) = e;
    elem(e,2) = e+1;
end
elem

elem = 100×2

1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11

K = sparse(N+1,N+1);
R = zeros(N+1,1);
k = 32;
L = 1;

for e = 1:N
    K(elem(e,1):elem(e,2), elem(e,1):elem(e,2)) = K(elem(e,1):elem(e,2), elem(e,1):elem(e,2)) + Kelem(e/N,N,L);
    R(elem(e,1):elem(e,2)) = R(elem(e,1):elem(e,2)) + Relem((elem(e,1)-1).*(L/N), (elem(e,2)-1).*(L/N),N,k,L);
end
function out = Kelem(out,N,L)
    Lelem = L/N;
    out = 1/2/Lelem*[integral(@(eta) out,-1,1,'ArrayValued',1) integral(@(eta) -out,-1,1,'ArrayValued',1); integral(@(eta) -out,-1,1,'ArrayValued',1) integral(@(eta) out,-1,1,'ArrayValued',1)];
end

function out = Relem(x1,x2,N,k,L)
    Lelem = L/N;
    phi1 = @(eta) 1/2*(1-eta);
    phi2 = @(eta) 1/2*(1+eta);
    out = [integral(@(eta) phi1(eta).*fun(eta,x1,x2,k,L), -1, 1).*(Lelem/2.0); integral(@(eta) phi2(eta).*fun(eta,x1,x2,k,L), -1, 1).*(Lelem/2.0)];
end

function out = fun(eta,x1,x2,k,L)
    out = -k^2*cos(2*pi*k*parametric(eta,x1,x2)/L);
end

function out = parametric(eta,x1,x2)
    out = (1-eta)/2.0*x1 + (1+eta)/2.0*x2;
end

function out = A(x)
    if x < .1
        A = @(x) 2 + 0*fun(eta,x1,x2,k,L);
    elseif x < .2
        A = @(x) 1.5 + 0*fun(eta,x1,x2,k,L);
    elseif x < .3
        A = @(x) 1.75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .4
        A = @(x) 1.5 + 0*fun(eta,x1,x2,k,L);
    elseif x < .5
        A = @(x) 3.75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .6
        A = @(x) .75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .7
        A = @(x) 1.25 + 0*fun(eta,x1,x2,k,L);
    elseif x < .8
        A = @(x) .75 + 0*fun(eta,x1,x2,k,L);
    elseif x < .9
        A = @(x) 2 + 0*fun(eta,x1,x2,k,L);
    elseif x <= 1
        A = @(x) 1 + 0*fun(eta,x1,x2,k,L);
    end
    out = A;
end

[VBBV]

since the output with function handle is named as out this must be given as input to function Kelem instead of x