I recently came, by chance, across this interesting post from one of the hpmuseum archives, originally posted by Valentin Albillo
**********************************************
Test 7  Solving a definite integral of an implicit function:
Find X in [0,1] such that
/X

 y(x).dx = 1/3

/0
where y(x) is an ultraradical function (a member of the family of elliptic functions) implicitly defined by:
5
y + y = x
using precisions of 1E3 and 1E6 for the integral.
HP71B code:
X=FNROOT(0,1,INTEGRAL(O,FVAR,1E3,FNROOT(0,1,FVAR^5+FVARIVAR))1/3)
this gives:
X = 0.854136725005 in 433 seconds (precision = 1E3)
= 0.854138746461 in 771 seconds (precision = 1E6)
************************************************
Obviously, this is very easy on the HP 71B, but it is also no problem for my HP 50G, with the following code
<< > x
<< 'Y^5+Y=x' 'Y' {0 1} ROOT >>
>>
'FY' STO
<<
<< 0 X 'FY(Z)' 'Z' \int 1 3 /  >>
'X' {0 1} ROOT
>>
EVAL
(and setting Number Format to Fix 3 or Fix 6 in the Mode Menu).
However, I run into problems with the HP Prime
Choosing 'FY' as a CAS program
(y) > BEGIN
RETURN(nSolve((t^5+t)=y , t));
END;
So far, so good: FY(x), x real, with nSolve gives me the one real solution, as I'm not interested in the four complex ones (which would pop up with, say, fsolve).
But now
nSolve(int(FY(z),z,0,x)1/3=0,x)
doesn't work, because FY can't take a symbolic argument (i.e., z). Or can it?
Any suggestions? There must be a simple way of doing this on the Prime, but my brain is stuck trying to translate this from HP50G syntax.
I don't believe there is a simple way to do that, since you don't have control on evaluation like in RPN.
Here is a not too complicated way in Xcas, therefore it should also work in CAS :
1/ first define the implict function
f(x):=fsolve(y^5+y=x,y=0..2)[0]
2/ then the integral g(t):=romberg('f(x)',x,0,t)
3/ Check that g(0.7) and g(1) have opposite sign.
4/ then solve using Newton method
h(t):=t(g(t)1/3)/f(t)
h(h(h(1.0))) and h(h(h(h(1.0)))) are almost identical.
Thank you very much, Bernard, this works!
May I ask three questions?
Step 1: why does this only work with fsolve and not with nSolve (here, the Prime hangs)?   by the way, the index on the Prime needs to be 1, not 0  0 is the convention in Xcas for the first element of a vector, I assume.
Step 2: why does this only work with romberg, but not with int?   also, by the way,  1/3 needs to be appended to the integral in step two, otherwise there will be no opposite signs in step 3, and then, in step 4, h(t): will be simply tg(t)/f(t).
Step 4: why does one have to resort to an explicit Newton construct to find the root, and cannot use fsolve again?
[P.S. Figured out why romberg works   because it returns the numerical result right away   so only two questions left!]
1/ nSolve was added in Xcas for TI compatibility, but the native command is fsolve, nSolve differs from fsolve in one important aspect: it does not autoquote argument, which is here important.
2/ I don't know why romberg works better than Gaussian quadratures... Perhaps this is again an evaluation problem
3/ Sign are indeed not opposite, I should have said check g(.7)1/3 and g(1)1/3, your definition of g is also correct.
4/ This is an evaluation problem, in algebraic notation you do not have enough control to make sure that the initial function f(x) will only be called with a numeric value for x. This can perhaps be improved by returning fsolve of the initial arguments.
Many thanks for the explanations!
Hi,
Why don't you use solve instead of fsolve?
Since y^5+y=x is a polynomial in y for each fixed x, solve(y^5+y=x,y) seems to be faster than fsolve.
Following up on my question, I wonder if you know why
fsolve((g(x))=1/3,x,0..1) returns no solutions.
Thanks
See above, you don't have enough control on evaluation to be able to use fsolve, because of the embedded fsolve and integration inside. They will generate errors, that are trapped, and after a finite number of errors, fsolve will return with [].
(06222014 06:31 AM)parisse Wrote: [ > ]See above, you don't have enough control on evaluation to be able to use fsolve, because of the embedded fsolve and integration inside. They will generate errors, that are trapped, and after a finite number of errors, fsolve will return with [].
thank you. I am still curious about the difference between solve and fsolve for polynomials. Is there a reference that you may provide?
The best reference is to look at the source code of giac in solve.cc, the in_fsolve function does the job :)
For polynomiallike equations, if you don't provide a guess or call bisection, proot is used. If the equation is not recognized to be polynomial then bisection is called.
(06232014 09:23 AM)parisse Wrote: [ > ]The best reference is to look at the source code of giac in solve.cc, the in_fsolve function does the job :)
For polynomiallike equations, if you don't provide a guess or call bisection, proot is used. If the equation is not recognized to be polynomial then bisection is called.
Thanks. It looks like the documentation for solve in Help is incomplete. It does not mention that you can place conditions on the unknowns:
solve(x^3+x+1=0) returns Warning! Algebraic extension not yet implemented for poly[1,0,1,1] after which it returns
the solutions in curly brackets {0.68232....complex soln's if complex is checked in CAS Settings}
solve(x^3+x+1=0, x=1..1) returns directly the real solution but within square brackets like so [0.68232...]. This usage of solve is not documented in Help (it seems to be just fsolve).
(06202014 02:40 PM)Helge Gabert Wrote: [ > ]I recently came, by chance, across this interesting post from one of the hpmuseum archives, originally posted by Valentin Albillo
**********************************************
Test 7  Solving a definite integral of an implicit function:
Find X in [0,1] such that
/X

 y(x).dx = 1/3

/0
where y(x) is an ultraradical function (a member of the family of elliptic functions) implicitly defined by:
5
y + y = x
using precisions of 1E3 and 1E6 for the integral.
Helge,
I apologize but I could not resist "hijacking" your post. It turns out that I am currently teaching Adv Calc.
and your question lead me to some good assignment for my students (some have the Prime). One problem (from Spivak's) was to show that there is a function f(x) such that (f(x))^5+(f(x))=x. Your question then leads to applications of the Fund Thm of Calc and to the Intermediate Value Thm. So, thanks!
BTW, it is really interesting to note how fast the Advanced Graphics of the HP Prime is for plotting the implicit function y^5+y=x compared to plotting F1(x):=solve(y^5+y=x,y) in the Function App
Alberto, you're very welcome   but Valentin Albillo is the one who originally came up with the post a few years back   I only dug it up and tried it on the Prime.
And yes, it is amazing how fast the implicit function plots in advanced graphing   it would be interesting to know what algorithm is used   it sure beats calling fsolve repeatedly! (It is also amazing to note how easy it is to program this problem on the 1980s era HP71B (with the optional Math module plugged in, I assume).
Very interesting   I will need to read both. Are you actually the author?
You should not be surprised that plotting an implicit function as such is faster than plotting it as an explicit function, where all calls to fsolve are independant. Hint: near a regular point you know the tangent, if you move a little along say x, you have a good starting point for the solver on the tangent. Singular points are harder to handle. This is the way plotimplicit works, it finds a discretization of the implicit curve, with all points of the discretization approx. solutions to the equation, once computed you can zoom, or compute approx area inside without recomputing. Another method if you want just a plot output would be to compute the sign pixel by pixel. It's probably possible to mix methods. Perhaps jte can summarize how the adv. graph works.
(06242014 05:26 AM)parisse Wrote: [ > ]You should not be surprised that plotting an implicit function as such is faster than plotting it as an explicit function, where all calls to fsolve are independant. Hint: near a regular point you know the tangent, if you move a little along say x, you have a good starting point for the solver on the tangent. Singular points are harder to handle. This is the way plotimplicit works, it finds a discretization of the implicit curve, with all points of the discretization approx. solutions to the equation, once computed you can zoom, or compute approx area inside without recomputing. Another method if you want just a plot output would be to compute the sign pixel by pixel. It's probably possible to mix methods. Perhaps jte can summarize how the adv. graph works.
Well, thanks. Did not know about plotimplcit. It is very interesting to see the two methods, plotimplicit and Jeff T advanced graphics work on this pitchfork (xy)^5+(xy)=(x+y)*(xy) when zooming.
I have just updated Xcas binaries, where you can now run
Code:
fsolve(int(fsolve(x^5+x=y,x,1e14),y,0,t)=1/3,t=1.0);
(fsolve returns itself if the equation has more than 1 variable)
(06202014 02:40 PM)Helge Gabert Wrote: [ > ]I recently came, by chance, across this interesting post from one of the hpmuseum archives, originally posted by Valentin Albillo
**********************************************
Test 7  Solving a definite integral of an implicit function:
Find X in [0,1] such that
/X

 y(x).dx = 1/3

/0
where y(x) is an ultraradical function (a member of the family of elliptic functions) implicitly defined by:
5
y + y = x
using precisions of 1E3 and 1E6 for the integral.
I took the liberty to implement it on a WP 34S. Since version 3.0 solve and integrate can be freely nested. The solution is relatively straight forward:
Code:
; y(x) using solver for label 00
LBL A
STO 00
# 000
# 001
SLV 00
RTN
; y^5 + y = x
LBL 00
x²
x²
INC X
×
RCL 0

RTN
; Integral  1/3
LBL B
PSE 0
# 000
RCL Y
⎰ A
# 003
1/x

RTN
# Solve label B with timing
LBL C
TICKS
STO 01
# 000
# 001
SLV B
TICKS
RCL 01
END
In FIX 03 mode, the result is 0.8540848... after 16.3 seconds.
In FIX 06 mode, the result is 0.8541387519... after 138.5 seconds.
The battery symbol flashes indicating that the button cells aren't fresh, probably degrading performance. It's still not bad. Interestingly, the time in FIX 06 mode is roughly ten times the time in FIX 03 mode. On the 71B higher accuracy just doubles the execution time. I guess this has to do with different integration and solving algorithms or at least a different handling of the requested accuracy.