Operands to the || and && operators must be convertible to logical scalar values in integrating product of functions using integral.

FiniteElement.m

classdef FiniteElement
    %FINITEELEMENT Summary of this class goes here
    
    properties
        N%number of nodes
        H% step in x direction
        domain % 0 to 1
        beta% hyperparameter
        
        %initial conditions
        alpha
        
        %analytic solution
        exact
        
        f%f(x, t)
        
        
        %Define boundary conditions
        g0  % alpha0
        g1  %alphaN
        g0dash  %alpha'0
        g1dash  %alpha'N
        
        %stiffness matrices
        A
        D
        
    end
    
    methods
        
        function obj = FiniteElement(N, beta)
            
            % Constructor initialize all variables
            %Parameters Initialization
            obj.N = N;
            obj.beta = beta;
            obj.H = 1 / N; % step in x direction
            obj.domain = linspace(0, 1, N + 1);
            obj.f = @(x, t) (1 - pi^2 * obj.beta ) * exp(t) * sin(pi * x);
            
            %Set initial conditions here
            obj.alpha = sin(pi * obj.domain(2:end - 1));
            
            %Set exact solution here
            obj.exact = @(t) exp(t) * sin(pi*obj.domain);
            %Define boundary conditions
            obj.g0 = @(t) 0;
            obj.g1 = @(t) 0;
            obj.g0dash = @(t) 0;
            obj.g1dash = @(t) 0;
            obj = obj.matrices_req();
        end
        
        function output = phi(obj, i, x)
            % basis functions see live script for better representation
            if i == 0
                
                if x >= 0 && x <= obj.H
                    output = -(x - 1 * obj.H) / obj.H;
                else
                    output = 0;
                end
                
            elseif i == obj.N
                
                if x >= (obj.N - 1) * obj.H  &&  x <= obj.N * obj.H
                    output = (x - (obj.N - 1) * obj.H) / obj.H;
                else
                    output = 0;
                end
                
            else
                
                if (x >= (i - 1) * obj.H) && (x <= i * obj.H)
                    output = (x - (i - 1) * obj.H) / obj.H;
                elseif x >= i * obj.H && x <= (i + 1) * obj.H
                    output = -(x - (i + 1) * obj.H) / obj.H;
                else
                    output = 0;
                end
                
            end
            
        end
        
        function output = phidash(obj, i, x)  % derivative
            
            %derivative of phi
            if i == 0
                
                if x >= 0 && x <= obj.H
                    output = -1 / obj.H;
                else
                    output = 0;
                end
                
            elseif i == obj.N
                
                if x >= (obj.N - 1) * obj.H && x <= obj.N * obj.H
                    output = 1 / obj.H;
                else
                    output = 0;
                end
                
            else
                
                if x >= (i - 1) * obj.H && x <= i * obj.H
                    output = 1 / obj.H;
                elseif x >= i * obj.H && x <= (i + 1) * obj.H
                    output = -1 / obj.H;
                else
                    output = 0;
                end
                
            end
            
        end
        
        function output = const(obj, i, t)       
            % constant in the matrix equation
            output = obj.g0dash(t) * integral(@(x) obj.phi(i, x) * obj.phi(0, x), 0, 1) + ...
                obj.g1dash(t) * integral(@(x) obj.phi(i, x) * obj.phi(obj.N, x), 0, 1) ...
                -obj.beta * obj.g0(t) * integral(@(x) obj.phidash(i, x) * obj.phidash(0, x), 0, 1) + ...
                -obj.beta * obj.g1(t) * integral(@(x) obj.phidash(i, x) * obj.phidash(obj.N, x), 0, 1)...
                +obj.beta*obj.phi(i, 1)*obj.phidash(obj.N, 1)*obj.g1(t)-obj.beta*obj.phi(i, 0)*obj.phidash(0, 0)*obj.g0(t) ;
        end
        
        function obj = matrices_req(obj)
            % stiffness matrices
            obj.A = zeros(obj.N - 1);
            obj.D = zeros(obj.N - 1);
            
            for i = 1:obj.N - 1
                
                for j = 1:obj.N - 1
                    
                    if ismember(abs(i - j), [0, 1])
                        obj.A(i, j) = integral(@(x) obj.phi(i, x) * obj.phi(j, x), 0, 1);
                        obj.D( i, j) = integral(@(x) obj.phidash(i, x) * obj.phidash(j, x), 0, 1);
                    end
                    
                end
                
            end
            
        end
        
        
        function derivative_ = dydx(obj, t, y)
            % derivative used in ode45
            C = zeros(obj.N - 1, 1);
            F = zeros(obj.N - 1, 1);
            
            for i = 1:obj.N - 1
                C(i) = obj.const(i, t);
                F(i) = integral(@(x) obj.phi(i, x) * obj.f(x, t), 0, 1);
                
            end
            
            derivative_ = obj.A \ (-(-obj.beta * obj.D * y + C) + F);
            
        end
        
        function output = exactsoln(obj, times)
            % give exactsolution the same structure as output from ode45
            
            output = zeros(length(times), length(obj.domain));
            
            for i = 1:length(times)
                output(i,:)  = obj.exact(times(i)) ;
            end
            
        end
        function [timestamps, ApproxTemp, ExactTemp] = temperature(obj, time_at)
            
            % main function to be called to get temperatures
            [timestamps, ApproxTemp] = ode45(@(t, y) obj.dydx(t, y), time_at, obj.alpha);
            ApproxTemp = horzcat( arrayfun(@(x) obj.g0(x),timestamps),...
                                ApproxTemp,...
                                   arrayfun(@(x) obj.g1(x), timestamps) );
            ExactTemp = obj.exactsoln(timestamps);
            
        end 
    end
    
end

run.m

obj = FiniteElement(15, -1);
[timestamp, A, E] = obj.temperature(linspace(.01, .1, 10)); 
x = obj.domain;
y = timestamp;
[xx, yy] = meshgrid(x, y);
subplot(2, 1, 1);
z = A;
heatmap(x, y, z);
colormap(flipud(hot));
title("Approximate Solution");
xlabel("The rod as x-axis");
ylabel("Time");



subplot(2, 1, 2);
z = E;
heatmap(x, y, z);
colormap(flipud(hot));
title("Exact Solution");
xlabel("The rod as x-axis");
ylabel("Time");

Error:


Error in FiniteElement/phi (line 76)
                if (x >= (i - 1) * obj.H) && (x <= i * obj.H)

Error in FiniteElement>@(x)obj.phi(i,x)*obj.phi(j,x) (line 140)
                        obj.A(i, j) = integral(@(x) obj.phi(i, x) * obj.phi(j, x), 0, 1);

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

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 88)
Q = integralCalc(fun,a,b,opstruct);

Error in FiniteElement/matrices_req (line 140)
                        obj.A(i, j) = integral(@(x) obj.phi(i, x) * obj.phi(j, x), 0, 1);

Error in FiniteElement (line 53)
            obj = obj.matrices_req();
 
>> ylabel("Time");

[Vishesh Mangla]

Problem is product of two functions.

Answer 1

[Ameer Hamza]

It happens if one of the operands is a vector. && only work when variables on both sides are scalar. Try your code after replacing && (double &) with & (single & operator).

[Vishesh Mangla]

That doesn't work. I had done discussions on the same answer with sorry I forgot the name of the MVP, but he has the icon with a "polar bear". He deleted his answer and I had given everything like errors and stuff and the proof. The problem is might in the internal implementation where instead of for loops some vectorization shortcut is used. I used https://www.mathworks.com/matlabcentral/fileexchange/28726-simpson-s-rule-integration in a prior code where @(x) obj.phi(i, x) * obj.phi(j, x) function was integrated as now and it too gave errors and also please don't expect me to do all of that work again. Please see the deleted answers of this question by other guys.

[Ameer Hamza]

Replacing && with & does resolve the mentioned error. But It gives a different error. It will be faster if you can create functions obj.phi such that it can accept a vector 'x' and return a vector. Right now, the intergal() function complains that is function cannot operate on a vector. As a quick-fix, you can specify that the integrand is arrayValued to the integral() function. Compare the following code with the current code and see the difference

classdef FiniteElement
    %FINITEELEMENT Summary of this class goes here
    
    properties
        N%number of nodes
        H% step in x direction
        domain % 0 to 1
        beta% hyperparameter
        
        %initial conditions
        alpha
        
        %analytic solution
        exact
        
        f%f(x, t)
        
        
        %Define boundary conditions
        g0  % alpha0
        g1  %alphaN
        g0dash  %alpha'0
        g1dash  %alpha'N
        
        %stiffness matrices
        A
        D
        
    end
    
    methods
        
        function obj = FiniteElement(N, beta)
            
            % Constructor initialize all variables
            %Parameters Initialization
            obj.N = N;
            obj.beta = beta;
            obj.H = 1 / N; % step in x direction
            obj.domain = linspace(0, 1, N + 1);
            obj.f = @(x, t) (1 - pi^2 * obj.beta ) * exp(t) * sin(pi * x);
            
            %Set initial conditions here
            obj.alpha = sin(pi * obj.domain(2:end - 1));
            
            %Set exact solution here
            obj.exact = @(t) exp(t) * sin(pi*obj.domain);
            %Define boundary conditions
            obj.g0 = @(t) 0;
            obj.g1 = @(t) 0;
            obj.g0dash = @(t) 0;
            obj.g1dash = @(t) 0;
            obj = obj.matrices_req();
        end
        
        function output = phi(obj, i, x)
            % basis functions see live script for better representation
            if i == 0
                
                if x >= 0 & x <= obj.H
                    output = -(x - 1 * obj.H) / obj.H;
                else
                    output = 0;
                end
                
            elseif i == obj.N
                
                if x >= (obj.N - 1) * obj.H  &  x <= obj.N * obj.H
                    output = (x - (obj.N - 1) * obj.H) / obj.H;
                else
                    output = 0;
                end
                
            else
                
                if (x >= (i - 1) * obj.H) & (x <= i * obj.H)
                    output = (x - (i - 1) * obj.H) / obj.H;
                elseif x >= i * obj.H & x <= (i + 1) * obj.H
                    output = -(x - (i + 1) * obj.H) / obj.H;
                else
                    output = 0;
                end
                
            end
            
        end
        
        function output = phidash(obj, i, x)  % derivative
            
            %derivative of phi
            if i == 0
                
                if x >= 0 & x <= obj.H
                    output = -1 / obj.H;
                else
                    output = 0;
                end
                
            elseif i == obj.N
                
                if x >= (obj.N - 1) * obj.H & x <= obj.N * obj.H
                    output = 1 / obj.H;
                else
                    output = 0;
                end
                
            else
                
                if x >= (i - 1) * obj.H & x <= i * obj.H
                    output = 1 / obj.H;
                elseif x >= i * obj.H & x <= (i + 1) * obj.H
                    output = -1 / obj.H;
                else
                    output = 0;
                end
                
            end
            
        end
        
        function output = const(obj, i, t)       
            % constant in the matrix equation
            output = obj.g0dash(t) * integral(@(x) obj.phi(i, x) * obj.phi(0, x), 0, 1, 'ArrayValued', 1) + ...
                obj.g1dash(t) * integral(@(x) obj.phi(i, x) * obj.phi(obj.N, x), 0, 1, 'ArrayValued', 1) ...
                -obj.beta * obj.g0(t) * integral(@(x) obj.phidash(i, x) * obj.phidash(0, x), 0, 1, 'ArrayValued', 1) + ...
                -obj.beta * obj.g1(t) * integral(@(x) obj.phidash(i, x) * obj.phidash(obj.N, x), 0, 1, 'ArrayValued', 1)...
                +obj.beta*obj.phi(i, 1)*obj.phidash(obj.N, 1)*obj.g1(t)-obj.beta*obj.phi(i, 0)*obj.phidash(0, 0)*obj.g0(t) ;
        end
        
        function obj = matrices_req(obj)
            % stiffness matrices
            obj.A = zeros(obj.N - 1);
            obj.D = zeros(obj.N - 1);
            
            for i = 1:obj.N - 1
                
                for j = 1:obj.N - 1
                    
                    if ismember(abs(i - j), [0, 1])
                        obj.A(i, j) = integral(@(x) obj.phi(i, x) * obj.phi(j, x), 0, 1, 'ArrayValued', 1);
                        obj.D( i, j) = integral(@(x) obj.phidash(i, x) * obj.phidash(j, x), 0, 1, 'ArrayValued', 1);
                    end
                    
                end
                
            end
            
        end
        
        
        function derivative_ = dydx(obj, t, y)
            % derivative used in ode45
            C = zeros(obj.N - 1, 1);
            F = zeros(obj.N - 1, 1);
            
            for i = 1:obj.N - 1
                C(i) = obj.const(i, t);
                F(i) = integral(@(x) obj.phi(i, x) * obj.f(x, t), 0, 1);
                
            end
            
            derivative_ = obj.A \ (-(-obj.beta * obj.D * y + C) + F);
            
        end
        
        function output = exactsoln(obj, times)
            % give exactsolution the same structure as output from ode45
            
            output = zeros(length(times), length(obj.domain));
            
            for i = 1:length(times)
                output(i,:)  = obj.exact(times(i)) ;
            end
            
        end
        function [timestamps, ApproxTemp, ExactTemp] = temperature(obj, time_at)
            
            % main function to be called to get temperatures
            [timestamps, ApproxTemp] = ode45(@(t, y) obj.dydx(t, y), time_at, obj.alpha);
            ApproxTemp = horzcat( arrayfun(@(x) obj.g0(x),timestamps),...
                                ApproxTemp,...
                                   arrayfun(@(x) obj.g1(x), timestamps) );
            ExactTemp = obj.exactsoln(timestamps);
            
        end 
    end
    
end

Running the script creates following graph in R2020a

[Vishesh Mangla]

Thanks @AmeerHamza , I 'm really surprised to see that this worked. I wasted almost 2 days in this error and finally used the simpsons38 method and submitted the assignment to my professor yesterday. I hope your answer will be useful to others.

[Vishesh Mangla]

Ahh, I missed you are having v2020a but I' was using v2018a. Maybe this might not work there.

[Ameer Hamza]

I think it should work in R2018b since these are very basic operators, there shouldn't be much difference. I got the same error after converting && to &. I have made some other changes too. See the integral() lines.

[Vishesh Mangla]

Yes, I did observe this

 'ArrayValued', 1

but I expected integral to evaluate values of the product of the phi's and sum them up like you do in trapazoidal or simpson's rule. I used the following code to submit my assignment https://www.mathworks.com/matlabcentral/fileexchange/72562-simpson-s-3-8-rule-composite and this gave me no errors.