how to convert Hypergeometric function generated in Maple to Matlab

4 views (last 30 days)
% Hey, I have an issue in converting below code in Matlab. I don't know how to write 'I1' in Matlab. if we run below code in Maple then plot will show decay with omega but that is not happening with Matlab due to error in I1 because of presence of "Hypergeom". Code is written in Maple 2018. Please help me in solving this issue. Thank You!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
restart; Digits:=500:
n:=4:
f1:=1: f2:=cos(alpha):
if n =0 then f3:=f1:
elif n =1 then f3:=f2:
else
for i from 3 by 1 to n+1 do f3:=expand(2*f2*cos(alpha)-f1);
f1:=f2: f2:=f3:
end do;
end if;
f3:
#I1_s:=simplify(int(x^k*sin(omega*x)/(1+x^2)^((n+k)/2+1),x=0..infinity)) assuming omega::real,omega>0,k::integer,k>0,n::integer,n>0;
#I1 := omega*hypergeom([1+(1/2)*k], [3/2, 1-(1/2)*n], (1/4)*omega^2)*GAMMA(1+(1/2)*k)*GAMMA((1/2)*n)/(2*GAMMA((1/2)*n+(1/2)*k+1))-omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*Pi/(2*sin((1/2)*n*Pi)*GAMMA(n+2));
#I1_c:=simplify(int(x^k*cos(omega*x)/(1+x^2)^((n+k)/2+1),x=0..infinity)) assuming omega::real,omega>0,k::integer,k>0,n::integer,n>0;
#I1 := (hypergeom([(1/2)*k+1/2], [1/2, 1/2-(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*k+1/2)*GAMMA((1/2)*n+1/2)*cos((1/2)*n*Pi)*GAMMA(n+2)-Pi*omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*n+(1/2)*k+1))/(2*GAMMA((1/2)*n+(1/2)*k+1)*cos((1/2)*n*Pi)*GAMMA(n+2));
## Depending on n=even or odd I1 is decided.
if irem(n,2)=0 then
I1 := (hypergeom([(1/2)*k+1/2], [1/2, 1/2-(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*k+1/2)*GAMMA((1/2)*n+1/2)*cos((1/2)*n*Pi)*GAMMA(n+2)-Pi*omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*n+(1/2)*k+1))/(2*GAMMA((1/2)*n+(1/2)*k+1)*cos((1/2)*n*Pi)*GAMMA(n+2));
elif irem(n,2)<>0 then
I1 := omega*hypergeom([1+(1/2)*k], [3/2, 1-(1/2)*n], (1/4)*omega^2)*GAMMA(1+(1/2)*k)*GAMMA((1/2)*n)/(2*GAMMA((1/2)*n+(1/2)*k+1))-omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*Pi/(2*sin((1/2)*n*Pi)*GAMMA(n+2));
end if;
h:=0:
for k from 0 by 1 to n do h:=h+(coeff(f3,cos(alpha),k)*I1):
end do:
s1:=evalf(h):
plot(s1,omega=0..100);return;

Answers (2)

Walter Roberson
Walter Roberson on 23 Dec 2021
Edited: Walter Roberson on 23 Dec 2021
hypergeom() is not the problem; it works the same in MATLAB. Just watch out for omega == 0 exactly
%digits(500);
syms alpha k omega
GAMMA = @gamma;
Pi = sym(pi);
n = 4;
f1 = 1; f2 = cos(alpha);
if n == 0
f3 = f1;
elsif n == 1
f3 = f2;
else
for i = 3 : n+1
f3 = expand(2*f2*cos(alpha)-f1);
f1 = f2; f2 = f3;
end
end
%f3:
%I1_s:=simplify(int(x^k*sin(omega*x)/(1+x^2)^((n+k)/2+1),x=0..infinity)) assuming omega::real,omega>0,k::integer,k>0,n::integer,n>0;
%I1 := omega*hypergeom([1+(1/2)*k], [3/2, 1-(1/2)*n], (1/4)*omega^2)*GAMMA(1+(1/2)*k)*GAMMA((1/2)*n)/(2*GAMMA((1/2)*n+(1/2)*k+1))-omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*Pi/(2*sin((1/2)*n*Pi)*GAMMA(n+2));
%I1_c:=simplify(int(x^k*cos(omega*x)/(1+x^2)^((n+k)/2+1),x=0..infinity)) assuming omega::real,omega>0,k::integer,k>0,n::integer,n>0;
%I1 := (hypergeom([(1/2)*k+1/2], [1/2, 1/2-(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*k+1/2)*GAMMA((1/2)*n+1/2)*cos((1/2)*n*Pi)*GAMMA(n+2)-Pi*omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*n+(1/2)*k+1))/(2*GAMMA((1/2)*n+(1/2)*k+1)*cos((1/2)*n*Pi)*GAMMA(n+2));
% Depending on n=even or odd I1 is decided.
if mod(n,2) == 0
I1 = (hypergeom([(1/2)*k+1/2], [1/2, 1/2-(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*k+1/2)*GAMMA((1/2)*n+1/2)*cos((1/2)*n*Pi)*GAMMA(n+2)-Pi*omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*n+(1/2)*k+1))/(2*GAMMA((1/2)*n+(1/2)*k+1)*cos((1/2)*n*Pi)*GAMMA(n+2));
elseif mod(n,2) ~= 0
I1 = omega*hypergeom([1+(1/2)*k], [3/2, 1-(1/2)*n], (1/4)*omega^2)*GAMMA(1+(1/2)*k)*GAMMA((1/2)*n)/(2*GAMMA((1/2)*n+(1/2)*k+1))-omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*Pi/(2*sin((1/2)*n*Pi)*GAMMA(n+2));
else
error('n is not a finite integer')
end
I1
I1 = 
f3
f3 = 
[C, pow] = coeffs(f3, cos(alpha), 'all')
C = 
pow = 
h = 0;
for K = 0 : min(n, length(C)-1)
h = h + C(end-K) * subs(I1, k, K);
end
h
h = 
s1 = vpa(h)
s1 = 
%fplot(s1, [0, 100]);
limit(s1, omega, 0)
ans = 
0.09817477042468103884879414609549
Omega = linspace(.1,100,250);
ds1 = double(subs(s1, omega, Omega));
plot(Omega, ds1)
plot(Omega(1:75), ds1(1:75))
  5 Comments
Walter Roberson
Walter Roberson on 24 Dec 2021
Edited: Walter Roberson on 24 Dec 2021
It turns out that the problems occur with some of the GAMMA and cos() and sin() calls; with double precision n values, those calls were operating at double precision resolution, which turned out not to be good enough.
digits(250);
syms alpha k omega
GAMMA = @gamma;
Pi = sym(pi);
n = 4;
N = sym(n);
f1 = 1; f2 = cos(alpha);
if n == 0
f3 = f1;
elseif n == 1
f3 = f2;
else
for i = 3 : n+1
f3 = expand(2*f2*cos(alpha)-f1);
f1 = f2; f2 = f3;
end
end
%f3:
%I1_s:=simplify(int(x^k*sin(omega*x)/(1+x^2)^((n+k)/2+1),x=0..infinity)) assuming omega::real,omega>0,k::integer,k>0,n::integer,n>0;
%I1 := omega*hypergeom([1+(1/2)*k], [3/2, 1-(1/2)*n], (1/4)*omega^2)*GAMMA(1+(1/2)*k)*GAMMA((1/2)*n)/(2*GAMMA((1/2)*n+(1/2)*k+1))-omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*Pi/(2*sin((1/2)*n*Pi)*GAMMA(n+2));
%I1_c:=simplify(int(x^k*cos(omega*x)/(1+x^2)^((n+k)/2+1),x=0..infinity)) assuming omega::real,omega>0,k::integer,k>0,n::integer,n>0;
%I1 := (hypergeom([(1/2)*k+1/2], [1/2, 1/2-(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*k+1/2)*GAMMA((1/2)*n+1/2)*cos((1/2)*n*Pi)*GAMMA(n+2)-Pi*omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*n+(1/2)*k+1))/(2*GAMMA((1/2)*n+(1/2)*k+1)*cos((1/2)*n*Pi)*GAMMA(n+2));
% Depending on n=even or odd I1 is decided.
if mod(n,2) == 0
I1 = (hypergeom([(1/2)*k+1/2], [1/2, 1/2-(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*k+1/2)*GAMMA((1/2)*N+1/2)*cospi((1/2)*N)*GAMMA(N+2)-Pi*omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*GAMMA((1/2)*n+(1/2)*k+1))/(2*GAMMA((1/2)*n+(1/2)*k+1)*cospi((1/2)*N)*GAMMA(N+2));
elseif mod(n,2) ~= 0
I1 = omega*hypergeom([1+(1/2)*k], [3/2, 1-(1/2)*n], (1/4)*omega^2)*GAMMA(1+(1/2)*k)*GAMMA((1/2)*N)/(2*GAMMA((1/2)*n+(1/2)*k+1))-omega^(1+n)*hypergeom([(1/2)*n+(1/2)*k+1], [1+(1/2)*n, 3/2+(1/2)*n], (1/4)*omega^2)*Pi/(2*sinpi((1/2)*N)*GAMMA(N+2));
else
error('n is not a finite integer')
end
[C, pow] = coeffs(f3, cos(alpha), 'all');
h = 0;
for K = 0 : min(n, length(C)-1)
part1 = C(end-K);
part2 = subs(I1, k, K);
h = h + part1 * part2;
%disp(K);
%disp(part1);
%disp(string(part2));
% disp(double(vpa(limit(part2, omega, 100),50)));
end
h
h = 
s1 = vpa(h)
s1 = 
Omega = linspace(.1,100,250);
ds1 = double(subs(s1, omega, Omega));
plot(Omega, ds1)

Sign in to comment.


gourav pandey
gourav pandey on 24 Dec 2021
hii, yes you have done lot of changes but still plot is incorrect.. its unbounded. it suppose to be bounded. i'm also getting same result as u have shared. i think there if some issue with hypergeom function in matlab

Categories

Find more on Event Functions in Help Center and File Exchange

Products


Release

R2012b

Community Treasure Hunt

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

Start Hunting!