Negative values while solving equation with gamma functions

8 次查看（过去 30 天）
Sreeja
Sreeja 2014-7-3
编辑： John D'Errico 2014-7-3
I am trying to solve the following equation which involves a gamma function. vpasolve(gamma(x+1.77)^2==gamma(x)*gamma(x+3.54)-gamma(x+1.77)^2,x), but the answer that I get is negative. Is there a way to solve this?

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

回答（1 个）

John D'Errico
John D'Errico 2014-7-3
编辑：John D'Errico 2014-7-3
Why do people feel the need to use vpasolve?
fzero will be faster, just not quite as accurate. But if your final goal is a double anyway, then why use a mack truck to take a pea to Boston? You can't store 40 digits of precision in a double.
Regardless, PLOT YOUR FUNCTION! Don't just throw up your hands and stop thinking!
ezplot(fun,[3,4])
Had you used a plot to look for a root, you would have found one. Then just give vpasolve a starting value.
vpasolve(gamma(x+1.77)^2==gamma(x)*gamma(x+3.54)-gamma(x+1.77)^2,x,4)
ans =
3.3472921715874174466992415057515
You could also have used a simple scheme to look for zero crossings programmatically, to then give you a bracket for the desired positive root.
  2 个评论
Sreeja
Sreeja 2014-7-3
编辑：Sreeja 2014-7-3
Thanks John.
Suppose negative values are obtained while running a for loop, then how can we solve it? For example,
y=[0.4:0.1:1]'; n=numel(y); syms x sym('result',size(y));
for j=1:n, result(j,1)= vpasolve(y(j)*y(j)*gamma(x+(2/1.0720))*gamma(x+(2/1.0720))==(gamma(x)*gamma(x+(4/1.0720))-(gamma(x+(2/1.0720))*gamma(x+(2/1.0720)))),x); end
Also which is the best effective way to solve this equation?
John D'Errico
John D'Errico 2014-7-3
编辑：John D'Errico 2014-7-3
As I suggested, generate a list of points of the function, then test for a sign change.
Yes, you DO have a function. Inside your loop, define a function handle as:
fun = @(x) -y(j)*y(j)*gamma(x+(2/1.0720)).*gamma(x+(2/1.0720)) + (gamma(x).*gamma(x+(4/1.0720))-(gamma(x+(2/1.0720)).*gamma(x+(2/1.0720))));
For a given value of y(i), this is a function of x. Evaluate at a list of points. Check to find the first consecutive pair of points such that
fun(x)*fun(x+dx) <= 0
Pass that as a bracket into fzero. Note that the test above can be computed in one line, with no need for a loop. Find will do most of the work.
Whats the problem?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by