how to convert Hypergeometric function generated in Maple to Matlab

% 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;

Answer 1

[Walter Roberson]

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))

[gourav pandey]

1stli i really appreciate your affort in converting the code and coming out with result.

but i have doubt in the result. actually in maple when we use "Digits:=500" then in the plot we see the decay as omega increase but here it show unbounded.

do we have something similar to :Digits" in Matlab. becuse without Digits command the maple also show the same result which you have shown.

Thank you.

[Walter Roberson]

I just took limit(h,omega,100) and copied it to Maple, and did a minor fixup on hypergeom syntax. The result evaluates the same in Maple as it does in MATLAB.

This tells us that Maple and MATLAB come up with different formulas for the integral.

[Walter Roberson]

A working version

digits(250);

syms alpha k omega
GAMMA = @gamma;
Pi = sym(pi);

n = 4;
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([(sym(1)/2)*k+sym(1)/2], [sym(1)/2, sym(1)/2-(sym(1)/2)*n], (sym(1)/4)*omega^2)*GAMMA((sym(1)/2)*k+sym(1)/2)*GAMMA((sym(1)/2)*n+sym(1)/2)*cos((sym(1)/2)*n*Pi)*GAMMA(n+2)-Pi*omega^(1+n)*hypergeom([(sym(1)/2)*n+(sym(1)/2)*k+1], [1+(sym(1)/2)*n, sym(3)/2+(sym(1)/2)*n], (sym(1)/4)*omega^2)*GAMMA((sym(1)/2)*n+(sym(1)/2)*k+1))/(2*GAMMA((sym(1)/2)*n+(sym(1)/2)*k+1)*cos((sym(1)/2)*n*Pi)*GAMMA(n+2));
elseif mod(n,2) ~= 0 
    I1 = omega*hypergeom([1+(sym(1)/2)*k], [sym(3)/2, 1-(sym(1)/2)*n], (sym(1)/4)*omega^2)*GAMMA(1+(sym(1)/2)*k)*GAMMA((sym(1)/2)*n)/(2*GAMMA((sym(1)/2)*n+(sym(1)/2)*k+1))-omega^(1+n)*hypergeom([(sym(1)/2)*n+(sym(1)/2)*k+1], [1+(sym(1)/2)*n, 3/2+(sym(1)/2)*n], (sym(1)/4)*omega^2)*Pi/(2*sin((sym(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)
  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

     0

1

(omega*pi*(6*cosh(omega) - (2*sinh(omega)*(omega^2 + 3))/omega))/32 - (3*pi*(omega*sinh(omega) - cosh(omega)*(omega^2/3 + 1)))/16

     0

     1

0

(90*pi^(1/2)*hypergeom(1, [-3/2, 1/2], omega^2/4) - (15*omega^5*pi^(3/2)*hypergeom([], 3, omega^2/4))/8)/(450*pi^(1/2))

  -5.8449e-05

     2

(pi*hypergeom(3/2, [-3/2, 1/2], omega^2/4))/32 - (omega^5*pi*hypergeom(4, [3, 7/2], omega^2/4))/240

   2.7866e-04

     3

0

(90*pi^(1/2)*hypergeom(2, [-3/2, 1/2], omega^2/4) - (105*omega^5*pi^(3/2)*hypergeom(9/2, [3, 7/2], omega^2/4))/16)/(1575*pi^(1/2))

  -6.0590e-04

     4

8

(3*pi*hypergeom(5/2, [-3/2, 1/2], omega^2/4))/256 - (omega^5*pi*hypergeom(5, [3, 7/2], omega^2/4))/240

   -0.0036

h

h = 

s1 = vpa(h)

s1 = 

        
%fplot(s1, [0, 100]);
limit(s1, omega, 0)

ans = 

0.09817477042468103870195760572748446513116154373047205690546701850961926269644403120712608829194115837444212770354007208272168264404845569474141748150401503492032137820943557878454873889446715484157443638840034298956042019151486400730583697391472505963

Omega = linspace(.1,100,250);
ds1 = double(subs(s1, omega, Omega));
plot(Omega, ds1)

%plot(Omega(1:75), ds1(1:75))

[Walter Roberson]

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)

[gourav pandey]

okay, thank you so much!!

Answer 2

[gourav pandey]

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