Solve with integrating an implicit function

[Helge Gabert]

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 ultra-radical function (a member of the family of elliptic functions) implicitly defined by:
5
y + y = x
using precisions of 1E-3 and 1E-6 for the integral.
HP-71B code:

X=FNROOT(0,1,INTEGRAL(O,FVAR,1E-3,FNROOT(0,1,FVAR^5+FVAR-IVAR))-1/3)
this gives:
X = 0.854136725005 in 433 seconds (precision = 1E-3)
= 0.854138746461 in 771 seconds (precision = 1E-6)

************************************************
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.

[parisse]

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.

[Helge Gabert]

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 t-g(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!]

[parisse]

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 auto-quote 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.

[Helge Gabert]

Many thanks for the explanations!

[Alberto Candel]

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.

[Alberto Candel]

Following up on my question, I wonder if you know why
fsolve((g(x))=1/3,x,0..1) returns no solutions.
Thanks

[parisse]

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 [].

[Alberto Candel]

(06-22-2014 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?

[parisse]

The best reference is to look at the source code of giac in solve.cc, the in_fsolve function does the job :-)
For polynomial-like 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.

[Alberto Candel]

(06-23-2014 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 polynomial-like 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).

[Alberto Candel]

(06-20-2014 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 ultra-radical function (a member of the family of elliptic functions) implicitly defined by:
5
y + y = x
using precisions of 1E-3 and 1E-6 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

[Helge Gabert]

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 HP-71B (with the optional Math module plugged in, I assume).

[jte]

(06-24-2014 01:01 AM)Helge Gabert Wrote:  ...implicit function plots in advanced graphing - - it would be interesting to know what algorithm is used...

If you'd like to learn more about how the Advanced Graphing app works, I suggest looking at

http://www.dgp.utoronto.ca/people/mooncake/msc.html

and / or

http://www.dgp.utoronto.ca/people/moonca...Tupper.pdf

[Helge Gabert]

Very interesting - - I will need to read both. Are you actually the author?

[parisse]

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.

[Alberto Candel]

(06-24-2014 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 (x-y)^5+(x-y)=(x+y)*(x-y) when zooming.

[parisse]

I have just updated Xcas binaries, where you can now run

Code:

fsolve(int(fsolve(x^5+x=y,x,1e-14),y,0,t)=1/3,t=1.0);

(fsolve returns itself if the equation has more than 1 variable)

[jte]

(06-24-2014 02:56 AM)jte Wrote:  http://www.dgp.utoronto.ca/people/mooncake/msc.html
http://www.dgp.utoronto.ca/people/moonca...Tupper.pdf

(06-24-2014 03:11 AM)Helge Gabert Wrote:  Very interesting - - I will need to read both.

(06-24-2014 03:11 AM)Helge Gabert Wrote:  Are you actually the author?

Yes.

[Marcus von Cube]

(06-20-2014 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 ultra-radical function (a member of the family of elliptic functions) implicitly defined by:
5
y + y = x
using precisions of 1E-3 and 1E-6 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.