how to convert Hypergeometric function generated in Maple to Matlab

3 ビュー (過去 30 日間)
gourav pandey
gourav pandey 2021 年 12 月 23 日
コメント済み: 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;

回答 (2 件)

Walter Roberson
Walter Roberson 2021 年 12 月 23 日
編集済み: Walter Roberson 2021 年 12 月 23 日
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 件のコメント
gourav pandey
gourav pandey 2021 年 12 月 24 日
okay, thank you so much!!

サインインしてコメントする。


gourav pandey
gourav pandey 2021 年 12 月 24 日
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

Community Treasure Hunt

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

Start Hunting!

Translated by