How to plot fsolve on a variable?
cadull
I require inverse solutions to the following function, but have difficulties
with both fsolve and solve.
> f := a -> -arcsin(0.303*sin(a)) + 13.8*a^2 - 24.5*a + 10.3;
> plot(f, 0.69..0.867); # plot domain of interest
fsolve gives correct answers, e.g.
> fsolve(f(a)=-0.5); # gives 0.745
solve gives a RootOf solution which when evaluated gives smaller values,
e.g.
> evalf(eval(solve(f(a)=t, a), t=-0.5)); # gives 0.0588
After passing through a few more functions, I use the results to plot a
graph. I cannot use fsolve because plot evaluates my function with a
variable. I cannot use solve because it isn't giving the desired answer.
Possible solutions:
1) manipulate solve to provide better solutions
2) force plot to call fsolve for each point evaluated
3) change the function to a form solve works with
(3) is not preferable since the form of f may change considerably,
however -0.3054*sin(a) +13.8*a^2 - 24.5*a + 10.3 is a close approximation
that works with solve.
(2) may have lower performance but could be the most reliable and speed is
not an issue for my application - is there a way to do this? I've searched
the documentation but haven't seen a method yet. I.e. is there a way to
ensure evaluation of a function (or procedure) is delayed until numerical
input can be provided? In my case, these values would come from plot.
(1) I've tried several things, such as placing assumptions on the variables,
but without success.
Any help appreciated, thanks. I'm using Maple 9.5.
Regards,
cadull
Peter Pein
or plot([f(a),a,a=1/2..1]);
--
Peter Pein
Berlin
Peter Pein
> Hi,
>
> I require inverse solutions to the following function, but have difficulties
> with both fsolve and solve.
>
>
>>f := a -> -arcsin(0.303*sin(a)) + 13.8*a^2 - 24.5*a + 10.3;
>>plot(f, 0.69..0.867); # plot domain of interest
> Any help appreciated, thanks. I'm using Maple 9.5.
>
> Regards,
> cadull
>
>
--
Peter Pein
Berlin
dev...@math.udel.edu
> I require inverse solutions to the following function, but have
difficulties
> with both fsolve and solve.
It is better to think of fsolve, solve, and RootOf as three distinct
solvers with fsolve and RootOf being the one that can give specific
solutions.
> > f := a -> -arcsin(0.303*sin(a)) + 13.8*a^2 - 24.5*a + 10.3;
> > plot(f, 0.69..0.867); # plot domain of interest
Note that you used the functional form of plot here: plot(f, ...)
rather than plot(f(a), a= ...). Why didn't you do that with the fsolve
has well? Had you done that, it would have worked.
I am going to define A:= 0.69; B:= 0.867; to make the rest of this
clearer. So the plot statement becomes
plot(f, A..B);
> fsolve gives correct answers, e.g.
> > fsolve(f(a)=-0.5); # gives 0.745
It is not the only correct answer, you just call it correct because it
happens to be the answer that you want. It is just random luck that
fsolve choses this answer from all the possibilities.
You can force fsolve to only give answers in a certain range:
f_inv:= x-> fsolve(f(a)=x, a, a= A..B);
plot(f_inv, f(A)..f(B));
> solve gives a RootOf solution which when evaluated gives smaller
values
Usually it does not do any good to try to restrict variable values
while using solve. Occasionally, something like
solve(f(a)=x, a) assuming a::RealRange(A,B);
will work, but it does not change the result in this case.
There is also a chance that solve will give a variety of analytic
solutions, one of which is the one you want, but it doesn't make any
difference in this case. It might help to set
_EnvAllSolutions:= true;
before calling solve, but it doesn't make any difference in this case.
On the other hand, range restrictions do work with RootOf, and RooOf
can be called directly, in lieu of solve.
f_inv:= x-> RootOf(f(a)=x, a, A..B);
plot(f_inv, f(A)..f(B));
Sometimes it is useful to have solve generate the RootOf's for you,
because they will be in a different form, perhaps they will evaluate
better. Indeed, in this case, the RootOf returned by solve is in a
different form than the original equation. The RootOf's thus returned
will not have the range information, but you can add it like this:
ans:= solve(f(a)=x, a);
subsindets(ans, specfunc(algebraic, RootOf), Z-> RootOf(op(Z), A..B));
> Possible solutions:
> 1) manipulate solve to provide better solutions
I gave some ideas above for how to do that, when it works. But it
doesn't work in general. This is not the fault of solve. There is
usually nothing analytic that can be done with a random or general
trancendental equation.
> 3) change the function to a form solve works with
It is very rare, in my experience, that that is useful.
> (2) may have lower performance but could be the most reliable and
speed is
> not an issue for my application - is there a way to do this? I've
searched
> the documentation but haven't seen a method yet. I.e. is there a way
to
> ensure evaluation of a function (or procedure) is delayed until
numerical
> input can be provided? In my case, these values would come from plot.
In general, you delay evaluation of any expression by using single
quotes:
plot('fsolve(...)', x= ...)
But I prefer not using the quotes if it can be avoided, and it can be
avoided in this case by using the functional form of plot.
There is a fourth alternative which is useful if you want to evaluate
the inverse rapidly and repeatedly: Generate a procedure which is
specialized to performing the Newton interation for this.
Dfa:= D(f);
finv:= proc(x)
# Newton interation for convex decreasing function
local an,r,oldDigits,newDigits,resn,res,Dfn;
if nargs=2 then oldDigits:= args[2] else oldDigits:= Digits fi;
Digits:= Digits+2;
newDigits:= oldDigits+2;
an:= .67; # Left endpoint
resn:= evalf(f(an)-x);
Dfn:= Dfa(an);
if resn < 0 then error "x must be greater than %1", f(an) fi;
do
r:= an-resn/Dfn;
res:= evalf(f(r)-x);
Dfn:= Dfa(r);
if fnormal(res, newDigits)=0 or fnormal(Dfn)=0 then
return evalf[oldDigits](r)
fi;
if r<an then
Digits:= Digits+ilog2(Digits);
return procname(x, oldDigits)
fi;
an,resn:= r,res
od
end proc;
Now
plot(finv, f(A)..f(B))
will work very quickly and accurately.
cadull
to acheive the desired objectives. I am now using RootOf with a range to
provide an inverse function along with the functional form of plot.
Previously I had thought the difference between plot(f, A..B) and plot(f(a),
a=A..B) was purely notational. I had not tried the functional form because
the function finally plotted was in two variables, one of which I passed a
constant as in plot(g(a, a0), a=A..B).
I appreciate the detailed response, I now have greater insight into how to
use these methods. In particular, the fourth alternative was interesting as
I tend to overlook such possibilities in a program that provides so much.
Regards,
cadull
<dev...@math.udel.edu> wrote in message
news:1111270030....@l41g2000cwc.googlegroups.com...
dev...@math.udel.edu
> Previously I had thought the difference between plot(f, A..B) and
plot(f(a),
> a=A..B) was purely notational. I had not tried the functional form
because
> the function finally plotted was in two variables, one of which I
passed a
> constant as in plot(g(a, a0), a=A..B).
In that case, you can use the functional form like this:
plot(a-> g(a,a0), A..B);