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 ExchangeLoops and Conditional Statements についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by