How can I plot the graph for this function?

f(x)=symsum(2^(-m*x)*(gamma(m*x+j-1+1)/gamma(m*x-1+1))*h(j/m),j.0.Inf)

where h(j/m)=(j/m)^2+1, m=100.

Please help me out.

[Dyuman Joshi]

Are you sure it's gamma and not square root?

[Gufran Malik]

@Dyuman Joshi. Yes sir, it's gamma function.

[Dyuman Joshi]

The code for the sum is as follows. Though, in this form, it takes to much time to compute.

If you analyse it on pen-paper you might get an reduced form which is easier/faster to compute.

syms x j
m=100;
2^(-m*x)*symsum(gamma(m*x+j)*((j/m)^2+1)/(gamma(m*x)*2^j*factorial(j)),j,0,Inf)

[Gufran Malik]

Thanks @Dyuman Joshi but this is giving error.

[Dyuman Joshi]

It worked on my pc without any errors, but it took an eternity to run :')

(i3-5th gen, 8 gb ddr3 ram)

Answer 1

[Torsten]

m = 100;
h = @(x) x.^2+1;
x = 0:0.01:10.0;
tol = 1e-16;
for i = 1:numel(x)
    Lm(i) = fun_Lm(m,h,x(i),tol);
end
plot(x,Lm)

function value = fun_Lm(m,h,x,tol)
value = 0.0;
error = Inf;
j = 0;
summand = 1.0;
while error > tol
    value = value + summand;
    error = abs(summand);
    summand = summand * (m*x+j) * 1/2 * 1/(j+1) * h((j+1)/m)/h(j/m);
    j = j + 1;
end
value = value/2^(m*x) ;
end

[Gufran Malik]

@Torsten Thank you so much. It worked for me and helped me a lot.

Answer 2

[Gufran Malik]

How can we plot this function? @Torsten Please help me with this

[Walter Roberson]

@Gufran Malik

You woiuld need a fixed value for m in order to plot that, unless you wanted to plot a surface with a range of m values.

It simplifies surprisingly far

syms h(x)
syms j m s x positive

h(x) = x^2 + 1

h(x) = 

inner(j) = int(h(s), s, j/m, (j+1)/m)

inner(j) = 

outer = symsum( gamma(m*x + j) / (gamma(m*x) * 2^j * factorial(j)) * inner(j), j, 0, inf)

outer = 

f(x) = outer / 2^(m*x)

f(x) = 

fs(x) = simplify(f(x))

fs(x) = 

[Torsten]

m = 100;
H = @(lb,ub) (1/3*ub^3+ub) - (1/3*lb^3+lb);     
x = 0:0.01:10.0;
tol = 1e-16;
Lm = zeros(size(x));
for i = 1:numel(x)
    Lm(i) = fun_Lm(m,H,x(i),tol);
end
Lm_Walter_Roberson = (9*m*x+1+3*m^2+3*(m*x).^2)/(3*m^3);
plot(x,[Lm;Lm_Walter_Roberson])

function value = fun_Lm(m,H,x,tol)
value = 0.0;
error = Inf;
j = 0;
summand = H(0/m,1/m);
while error > tol
    value = value + summand;
    error = abs(summand);
    summand = summand * (m*x+j) * 1/2 * 1/(j+1) * H((j+1)/m,(j+2)/m)/H(j/m,(j+1)/m);   
    j = j + 1;
end
value = value/2^(m*x) ;
end

[Gufran Malik]

@Torsten Thank you so much.

[Gufran Malik]

@Walter Roberson Thank you so much.