can I find a solution for this integral?

r=sqrt(x.^2+y.^2);
theta=atan2(y,x);
ab1=1;xi=2;n=10;
h=@(t) (t.*(exp(-xi.*(r.*cos(theta)+ab1))-exp(-t.*(r.*cos(theta)+ab1))).*((-1).^(n-1).*exp(-ab1.*t)./factorial( n )).*t.*besselj(1,t.*r.*sin(theta)));
integral(h, 0, inf)

7 comentarios

Shreen El-Sapa
Shreen El-Sapa el 3 de Ag. de 2021
when i run, it gives
Arrays have incompatible sizes for this operation.
Error in
test>@(t1)((t1-sqrt(t1.^2+k.^2)).^(-1).*(t1.*(exp(-sqrt(t1.^2+k.^2).*(r.*cos(t)+ab1))-exp(-t1.*(r.*cos(t)+ab1))).*(((-1).^(i-1).*t1.^(i-1).*exp(-ab1.*t1)./factorial(i))+((-1).^(i).*sqrt(pi.*k./2/t1.^2).*exp(-ab1.*sqrt(t1.^2+k.^2)).*gegenbauerC(i,-1./2,sqrt(t1.^2+k.^2)./k))).*t1.*besselj(1,t1.*r.*sin(t))))
(line 78)
h=@(t1)((t1-sqrt(t1.^2+k.^2)).^(-1).*
(t1.*(exp(-sqrt(t1.^2+k.^2).*(r.*cos(t)+ab1))-exp(-t1.*(r.*cos(t)+ab1))).*(((-1).^(i-1).*t1.^(i-1).*exp(-ab1.*t1)./factorial( i
))+((-1).^(i).*sqrt(pi.*k./2/t1.^2).*exp(-ab1.*sqrt(t1.^2+k.^2)).*gegenbauerC(i,-1./2,
sqrt(t1.^2+k.^2)./k))).*t1.*besselj(1,t1.*r.*sin(t))));
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 87)
Q = integralCalc(fun,a,b,opstruct);
Error in test (line 81)
h11=integral(h, 0, 10); %inf
Related documentation
>>
Walter Roberson
Walter Roberson el 3 de Ag. de 2021
What is size(x), size(y), size(n) ?
Shreen El-Sapa
Shreen El-Sapa el 3 de Ag. de 2021
Editada: Walter Roberson el 3 de Ag. de 2021
A=[ -0.00767
0.00696
-0.00031
0.00005
-0.00001
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000];
B=[ 0.12000
0.00209
-0.00019
0.00001
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
0.00000
.00000E+00
.00000E+00
.00000E+00
.00000E+00
0.00000
0.00000
0.00000
0.00000
0.00000];
a = 1 ; %RADIUS
L=.4;
c =-a/L;
b =a/L;
m =a*200; % NUMBER OF INTERVALS
[x,y]=meshgrid((c:(b-c)/m:b),(c:(b-c)/m:b)');
[I J]=find(sqrt(x.^2+y.^2)<(a-.1));
if ~isempty(I)
x(I,J) = 0;
y(I,J) = 0;
end
r=sqrt(x.^2+y.^2);
t=atan2(y,x);
warning off
psi1=0;
chi=1;alpha=1;ab1=1;
k=sqrt(chi.^2+alpha.^2); %xi=sqrt(t.^2+k.^2);
for i=2:20
h=@(t1)((t1-sqrt(t1.^2+k.^2)).^(-1).* (t1.*(exp(-sqrt(t1.^2+k.^2).*(r.*cos(t)+ab1))-exp(-t1.*(r.*cos(t)+ab1))).*(((-1).^(i-1).*t1.^(i-1).*exp(-ab1.*t1)./factorial( i ))+((-1).^(i).*sqrt(pi.*k./2/t1.^2).*exp(-ab1.*sqrt(t1.^2+k.^2)).*gegenbauerC(i,-1./2, sqrt(t1.^2+k.^2)./k))).*t1.*besselj(1,t1.*r.*sin(t))));
v=@(t1)((t1-sqrt(t1.^2+k.^2)).^(-1).* ((-t1.*exp(-sqrt(t1.^2+k.^2).*(r.*cos(t)+ab1))+sqrt(t1.^2+k.^2).*exp(-t1.*(r.*cos(t)+ab1))).*(((-1).^(i).*t1.^(i-1).*exp(-ab1.*t1)./factorial( i ))+((-1).^(i-1).*sqrt(pi.*k./2/(t1.^2+k.^2)).*exp(-ab1.*sqrt(t1.^2+k.^2)).*gegenbauerC(i,-1./2, sqrt(t.^2+k.^2)./k))).*t1.*besselj(1,t1.*r.*sin(t))));
h11=integral(h, 0, 10); %inf
h22=integral(v, 0, 10);%inf
%h11=0;h22=0;
Ai=A(i-1);Bi=B(i-1);
psi1=psi1+(Ai.*(r.^(-i+1)+h11)+(r.^(1./2).*besselk(i-1./2,r.*alpha)+h22).*Bi).*gegenbauerC(i,-1./2, cos(t));
end
[DH1,hh2]=contour(x,y,psi,10,'-k');
%clabel(CH,h1,'FontSize',8);clabel(DH,h2,'FontSize',8)
hold on
m1=100;
r1=ones(1,m1+1)*a;
th=(0:2*pi/m1:2*pi);
set(polar(th,r1,'-k'),'LineWidth',1.1);
title('$\kappa=0$','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')
ylabel({'$\eta=1\quad$'},'Interpreter','latex','FontSize',10,'rot',360,'FontName','Times New Roman','FontWeight','Normal');
%title('Happel$^\prime$s model','Interpreter','latex','FontSize',10,'FontName','Times New Roman','FontWeight','Normal')
%axis square
axis off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
format long g
n = randi([5 10])
n =
10
x = randn(1,5)
x = 1×5
-1.81858516726184 -0.416052891266473 1.66042693477506 0.5746136191372 -1.46363015622315
y = randn(1,5)
y = 1×5
-2.10212774197359 1.3805995039459 0.291512314848963 -0.799668147618347 -1.30627881007708
syms t real
r = sqrt(x.^2+y.^2);
theta = atan2(y,x);
ab1=1; xi=2; n=10;
h=@(t) (t.*(exp(-xi.*(r.*cos(theta)+ab1))-exp(-t.*(r.*cos(theta)+ab1))).*((-1).^(n-1).*exp(-ab1.*t)./factorial( n )).*t.*besselj(1,t.*r.*sin(theta)));
hint = integral(h, 0, inf, 'ArrayValued', true)
Warning: Inf or NaN value encountered.
hint = 1×5
NaN 1.9484262526967e-08 3.60543882192953e-10 -3.713705199688e-09 NaN
Walter Roberson
Walter Roberson el 3 de Ag. de 2021
The NaN appears to show up when x and y are both negative.
Shreen El-Sapa
Shreen El-Sapa el 3 de Ag. de 2021
i will use this notation and recalculate it.
thanks so much
Shreen El-Sapa
Shreen El-Sapa el 4 de Ag. de 2021
in this case it excuted but i wanted to use the values in the lines 17 and 18 it failed
syms t real
a = 1 ; %RADIUS
L=.4;
c =-a/L;
b =a/L;
m =a*200; % NUMBER OF INTERVALS
[x,y]=meshgrid([c:(b-c)/m:b],[c:(b-c)/m:b]');
[I J]=find(sqrt(x.^2+y.^2)<(a-.1));
if ~isempty(I)
x(I,J) = 0;
y(I,J) = 0;
end
r = sqrt(x.^2+y.^2);
theta = atan2(y,x);
ab1=1; xi=2; n=10; k=2;
% I want to use the following green lines instead of above values of ab1=1; xi=2; n=10; k=2;
%chi=1;alpha=1;ab1=1;
%k=sqrt(chi.^2+alpha.^2);xi=sqrt(t.^2+k.^2);
h=@(t)(t-xi).^(-1).* (t.*(exp(-xi.*(r.*cos(theta)+ab1))-exp(-t.*(r.*cos(theta)+ab1))).*(((-1).^(n-1).*t.^(n-1).*exp(-ab1.*t)./factorial( n ))+((-1).^(n).*sqrt(pi.*k./2./t.^2).*exp(-ab1.*xi).*gegenbauerC(n,-1./2, xi./k))).*t.*besselj(1,t.*r.*sin(theta)));
hint = integral(h, 0, inf, 'ArrayValued', true)

Iniciar sesión para comentar.

 Respuesta aceptada

Dave B
Dave B el 3 de Ag. de 2021
When I run this code with scalar x,y it works fine:
x=2;y=2;
r=sqrt(x.^2+y.^2);
theta=atan2(y,x);
ab1=1;xi=2;n=10;
h=@(t) (t.*(exp(-xi.*(r.*cos(theta)+ab1))-exp(-t.*(r.*cos(theta)+ab1))).*((-1).^(n-1).*exp(-ab1.*t)./factorial( n )).*t.*besselj(1,t.*r.*sin(theta)));
integral(h, 0, inf)
When I run this code with vector x and y it fails:
x=[1 2];y=[1 2];
r=sqrt(x.^2+y.^2);
theta=atan2(y,x);
ab1=1;xi=2;n=10;
h=@(t) (t.*(exp(-xi.*(r.*cos(theta)+ab1))-exp(-t.*(r.*cos(theta)+ab1))).*((-1).^(n-1).*exp(-ab1.*t)./factorial( n )).*t.*besselj(1,t.*r.*sin(theta)));
integral(h, 0, inf)
"For scalar-valued problems, the function y = fun(x) must accept a vector argument, x, and return a vector result, y. This generally means that fun must use array operators instead of matrix operators. For example, use .* (times) rather than * (mtimes). If you set the 'ArrayValued' option to true, then fun must accept a scalar and return an array of fixed size."
Here, your function accepts a scalar argument (t) and returns an array:
size(h(0)) == size(r); % for all scalar h, for all size(r) (at least that I tested)
"Set this flag to true or 1 to indicate that fun is a function that accepts a scalar input and returns a vector, matrix, or N-D array output."
Thus, if x and y are arrays, do you want the following?
integral(h, 0, inf, 'ArrayValued', true)

2 comentarios

Shreen El-Sapa
Shreen El-Sapa el 3 de Ag. de 2021
thanks so much
Shreen El-Sapa
Shreen El-Sapa el 4 de Ag. de 2021
in this case it excuted but i wanted to use the values in the lines 17 and 18 it failed
syms t real
a = 1 ; %RADIUS
L=.4;
c =-a/L;
b =a/L;
m =a*200; % NUMBER OF INTERVALS
[x,y]=meshgrid([c:(b-c)/m:b],[c:(b-c)/m:b]');
[I J]=find(sqrt(x.^2+y.^2)<(a-.1));
if ~isempty(I)
x(I,J) = 0;
y(I,J) = 0;
end
r = sqrt(x.^2+y.^2);
theta = atan2(y,x);
ab1=1; xi=2; n=10; k=2;
% I want to use the following green lines instead of above values of ab1=1; xi=2; n=10; k=2;
%chi=1;alpha=1;ab1=1;
%k=sqrt(chi.^2+alpha.^2);xi=sqrt(t.^2+k.^2);
h=@(t)(t-xi).^(-1).* (t.*(exp(-xi.*(r.*cos(theta)+ab1))-exp(-t.*(r.*cos(theta)+ab1))).*(((-1).^(n-1).*t.^(n-1).*exp(-ab1.*t)./factorial( n ))+((-1).^(n).*sqrt(pi.*k./2./t.^2).*exp(-ab1.*xi).*gegenbauerC(n,-1./2, xi./k))).*t.*besselj(1,t.*r.*sin(theta)));
hint = integral(h, 0, inf, 'ArrayValued', true)

Iniciar sesión para comentar.

Más respuestas (0)

Preguntada:

el 3 de Ag. de 2021

Comentada:

el 4 de Ag. de 2021

Community Treasure Hunt

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

Start Hunting!

Translated by