i solved the problem by dividing it up in 4 integrals, since the N vector isnt big (2×1), so I calculated all 4 elements of the 2×2 matrix and then put them back together into a matrix
quadgk problems
2 views (last 30 days)
Show older comments
The complete problem is very complex, and as in my last question regarding the same project I skipped the introduction, and I will do so again here.
The specific problem is (this time) an Array (2 cells) of function handles, and I am getting the same error every time (regardless of the changes I make in it (* -> .*, A'*B -> (A*B)' and so on...). So here is the code:
z = sym ('z','real');
zd = sym ('zd','real');
%for demonstration the values of M,w,mi,eps... are simplified
M=4;mi=1;eps0=1;w=100;a=1;b=1;h=10;
deltaz=h/M;
%functions g, R & k
k = w * sqrt(mi * eps0);
R = @(z,zd) sqrt((z-zd).*(z-zd) + a.*a);
g = @(z,zd) (exp(-1i.*k.*R(z,zd))./R(z,zd));
for r = 1:1:M
for s = 1:1:M
%N for index s
Ni1 = @(z) ((s + 1)*deltaz - z)./deltaz;
Ni2 = @(z) (z - s*deltaz)./deltaz;
Ni = {Ni1 Ni2}';
%N for index j
Nj1 = @(z) ((r + 1)*deltaz - z)./deltaz;
Nj2 = @(z) (z - r*deltaz)./deltaz;
Nj = {Nj1 Nj2}';
%N'
Npi1 = @(zd) ((s + 1)*deltaz - zd)./deltaz;
Npi2 = @(zd) (zd - s*deltaz)./deltaz;
Npi = {Npi1 Npi2}';
%D = d(N) / d(z)
D1 = -1./deltaz;
D2 = 1./deltaz;
D = cell2mat({D1 D2})';
%D'
Dp1 = D1;
Dp2 = D2;
Dp = cell2mat({Dp1 Dp2})';
%FIRST PROBLEM
funkcija1 = @(z,zd)(Dp*g(z,zd))';
funkcija2 = @(z)D*quadgk(@(zd)funkcija1(z,zd),s*deltaz,(s+1)*deltaz);
IntPrvi = quadgk(funkcija2(z),r*deltaz, (r+1)*deltaz);
%SECOND PROBLEM
funk3 = @(z,zd) (Npi(zd)*g(z,zd))';
funk4 = @(z)Nj(z)*int(funk3,zd,s*deltaz, (s+1)*deltaz);
IntDrugi = k.*k.*quadgk(@(z)funk4(z), r*deltaz, (r+1)*deltaz);
%here comes a blok of code that uses IntPrvi & IntDrugi
end
end
for the first integral the problem is matrix sizes and return value type; TRACE:
_??? Error using ==> quadgk>finalInputChecks at 476 Supported classes are 'double' and 'single'.
Error in ==> quadgk>evalFun at 358 finalInputChecks(x,fx);
Error in ==> quadgk>f1 at 375 [y,too_close] = evalFun(tt);
Error in ==> quadgk>vadapt at 269 [fx,too_close] = f(x);
Error in ==> quadgk at 208 [q,errbnd] = vadapt(@f1,interval);
Error in ==> struja>@(z)D*quadgk(@(zd)funkcija1(z,zd),s*deltaz,(s+1)*deltaz) at 94 funkcija2 = @(z)D*quadgk(@(zd)funkcija1(z,zd),s*deltaz,(s+1)*deltaz); %D*int(funk1)
Error in ==> struja at 95 IntPrvi = quadgk(funkcija2(z),r*deltaz, (r+1)*deltaz); _ ??? Error using ==> quadgk>finalInputChecks at 476 Supported classes are 'double' and 'single'.
Error in ==> quadgk>evalFun at 358 finalInputChecks(x,fx);
Error in ==> quadgk>f1 at 375 [y,too_close] = evalFun(tt);
Error in ==> quadgk>vadapt at 269 [fx,too_close] = f(x);
Error in ==> quadgk at 208 [q,errbnd] = vadapt(@f1,interval);
Error in ==> struja>@(z)D*quadgk(@(zd)funkcija1(z,zd),s*deltaz,(s+1)*deltaz) at 94 funkcija2 = @(z)D*quadgk(@(zd)funkcija1(z,zd),s*deltaz,(s+1)*deltaz); %D*int(funk1)
Error in ==> struja at 95 IntPrvi = quadgk(funkcija2(z),r*deltaz, (r+1)*deltaz);_
for the second integral the problem is how to call the Function handle Array Ni/Nj with the variable z, and multiply it (matrix multipl.) with g(z,zd) or with the result of quadprog for funk4. TRACE: _ ??? Subscript indices must either be real positive integers or logicals.
Error in ==> struja>@(z)Nj(z)*int(funk3,zd,s*deltaz,(s+1)*deltaz)
Error in ==> struja>@(z)funk4(z) at 100 IntDrugi = k.*k.*quadgk(@(z)funk4(z), r*deltaz, (r+1)*deltaz);
Error in ==> quadgk>evalFun at 357 fx = FUN(x);
Error in ==> quadgk>f1 at 375 [y,too_close] = evalFun(tt);
Error in ==> quadgk>vadapt at 269 [fx,too_close] = f(x);
Error in ==> quadgk at 208 [q,errbnd] = vadapt(@f1,interval);
Error in ==> struja at 100 IntDrugi = k.*k.*quadgk(@(z)funk4(z), r*deltaz, (r+1)*deltaz); _
I would be VERY gratefull for any insight why these errors are popping up, or any suggestion for a solution
0 Comments
Accepted Answer
More Answers (2)
Walter Roberson
on 12 Jul 2011
For your second problem, your line
funk4 = @(z)Nj(z)*int(funk3,zd,s*deltaz, (s+1)*deltaz);
is one of the few places that does not bind both z and zd as the dummy argument names for the anonymous function. The zd referred to will thus be the sym() that you created with that name.
The int() of a function is a symbolic operation and will always return a symbolic value, so funk4 will return a symbolic value, which quadgk() will not be able to deal with. If you are sure that the integral will be fully resolved to a numeric value, you could wrap the int() in double() to get the MATLAB floating point equivalent of the symbolic number. If you are expecting the int() to return a symbolic function, you may need to use matlabFunction() to convert that symbolic expression in to an anonymous function that quadgk() could integrate over.
I have not worked through your first problem, but your line
IntPrvi = quadgk(funkcija2(z),r*deltaz, (r+1)*deltaz);
is going to be working with symbolic expressions it appears to me.
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!