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

1 view (last 30 days)
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");

Accepted Answer

Ameer Hamza
Ameer Hamza on 19 Jun 2020
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).
  6 Comments
Ameer Hamza
Ameer Hamza on 21 Jun 2020
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
Vishesh Mangla on 21 Jun 2020
Edited: Vishesh Mangla on 21 Jun 2020
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Debugging and Analysis in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!