how to convert Hypergeometric function generated in Maple to Matlab

Question

gourav pandey 2021-12-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1616390-how-to-convert-hypergeometric-function-generated-in-maple-to-matlab

评论： gourav pandey 2021-12-24

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2021-12-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1616390-how-to-convert-hypergeometric-function-generated-in-maple-to-matlab#answer_860700

编辑：Walter Roberson 2021-12-23

在 MATLAB Online 中打开

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

%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 个评论
显示 3更早的评论隐藏 3更早的评论

Walter Roberson 2021-12-24

编辑：Walter Roberson 2021-12-24

在 MATLAB Online 中打开

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

%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 2021-12-24

编辑：Walter Roberson 2021-12-24

在 MATLAB Online 中打开

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

%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 2021-12-24

okay, thank you so much!!

请先登录，再进行评论。

Answer 2

gourav pandey 2021-12-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1616390-how-to-convert-hypergeometric-function-generated-in-maple-to-matlab#answer_861360

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

how to convert Hypergeometric function generated in Maple to Matlab

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

5 个评论
显示 3更早的评论隐藏 3更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

how to convert Hypergeometric function generated in Maple to Matlab

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

5 个评论 显示 3更早的评论隐藏 3更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论