Problem using quad2d with vector

1 view (last 30 days)
Scott Lines
Scott Lines on 26 Oct 2016
Answered: Torsten on 26 Oct 2016
Using below code with the variable t doesn't work, matrix dimensions don't agree, if the variable is replaced with a single number it does work (though that's not particularly useful), I've tried replacing it with a loop to no avail.
What am i missing?
lm = 1.75;
Cm = 3.2E6;
fe = 1380;
H = 13.5;
% time and space
y = 0.0;
ro = 0.1;
x = ro;
ts = [0.1 97/24];
Mt = 100;
% t = logspace(log10(ts(1)*86400),log10(ts(2)*86400),Mt); % [s]
t = linspace(ts(1)*86400, ts(2)*86400, Mt); % [s]
% calculated parameters
QL = fe/H;
z = H/2;
Dt = lm/Cm;
r = sqrt(x.^2+y.^2);
T = QL*(1/4/pi/lm)*...
(quad2d(@(f,ze) erfc(sqrt(r^2+ro^2-2*r*ro*cos(f)+(z-ze).^2)...
/2/sqrt(Dt*t))/pi./sqrt(r^2+ro^2-2*r*ro*cos(f)+(z-ze).^2),0,pi,0,H)-...
quad2d(@(f,ze) erfc(sqrt(r^2+ro^2-2*r*ro*cos(f)+(z-ze).^2)...
/2/sqrt(Dt*t))/pi./sqrt(r^2+ro^2-2*r*ro*cos(f)+(z-ze).^2),0,pi,-H,0));

Accepted Answer

KSSV
KSSV on 26 Oct 2016
You have to take care of element wise operations in that big expression T. As T is very big and clumsy, I have modified it to a loop, as below.
clc; clear all ;
lm = 1.75;
Cm = 3.2E6;
fe = 1380;
H = 13.5;
% time and space
y = 0.0;
ro = 0.1;
x = ro;
ts = [0.1 97/24];
Mt = 100;
% t = logspace(log10(ts(1)*86400),log10(ts(2)*86400),Mt); % [s]
tt = linspace(ts(1)*86400, ts(2)*86400, Mt); % [s]
% calculated parameters
QL = fe/H;
z = H/2;
Dt = lm/Cm;
r = sqrt(x.^2+y.^2);
T = zeros(size(tt)) ;
for i = 1:length(tt)
t = tt(i) ;
T(i) = QL*(1/4/pi/lm)*...
(quad2d(@(f,ze) erfc(sqrt(r^2+ro^2-2*r*ro*cos(f)+(z-ze).^2)...
/2/sqrt(Dt*t))/pi./sqrt(r^2+ro^2-2*r*ro*cos(f)+(z-ze).^2),0,pi,0,H)-...
quad2d(@(f,ze) erfc(sqrt(r^2+ro^2-2*r*ro*cos(f)+(z-ze).^2)...
/2/sqrt(Dt*t))/pi./sqrt(r^2+ro^2-2*r*ro*cos(f)+(z-ze).^2),0,pi,-H,0));
end
Though it is working, it is throwing some warning. I suggest you to check the expression for T.
  1 Comment
Scott Lines
Scott Lines on 26 Oct 2016
Thank you, I think I will need to go over the expression again but that's solved the first problem.

Sign in to comment.

More Answers (1)

Torsten
Torsten on 26 Oct 2016
There is no other way than to loop over the elements of t and call quad2d for each element separately.
MATLAB's "integral" has the option "ArrayValued=True", but as far as I know, integral2 or quad2d haven't.
Best wishes
Torsten.

Tags

Community Treasure Hunt

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

Start Hunting!