Hello everyone,

Here is what confuses me: There is a function with the symbol "I" in it, and Maple can still plot it in the 2-d Euclidean space. And here are the details:

The function I want to find the inverse of is f := x->piecewise(x <= 2, 13-8*x+4*x^2-x^3, x <= 3, -39.13953960+64.18605280*x-29.08139810*x^2+4.011628300*x^3, x <= 9, 8.013637628-3.472493940*x+.3065739839*x^2-0.1316006154e-1*x^3)

I wrote a procedure (with great help from Robert Israel of mapleprimes.com) to invert this function and the result I get is f_inverse:=y->Invert(2, 13-8*x+4*x^2-x^3, 3, -39.13953960+64.18605280*x-29.08139810*x^2+4.011628300*x^3, 9, 8.013637628-3.472493940*x+.3065739839*x^2-0.1316006154e-1*x^3):

where the procedure Invert() is pasted below.

When I plot f(x) and f_inverse(y), only in the first quadrant, they really look like inverses. But if you take time to look at the equation defining f_inverse(y), you will see I's in it.

First of all, I do not understand how Invert() returns an expression with complex variables in it. I thought I controlled for it.

Second, how can maple plot it with the command plot? Later, I need to assign diff(y*f_inverse(y),y) to another function and plot that as well. Same problem there...

You can paste the code below to a worksheet and execute.

Any comments clarifying these questions will be greatly appreciated.

Thank you all in advance.

Ozgur Inal

Invert := proc()
#enter the piecewise function to be inverted, say, piecewise(p<=1, p^2, p<=4, p^3) as InverseFunction(1,p^2,4,p^3).
local i, j, n, f, lx, ly, L, P, finverse, RD_q, MR;
#lp and lq are the coordinates of the points of the original function at the horizontal and vertical axis, respectively.
n := nargs/2;                                  
lx[1] := 0;
ly[1] := 0;
#initialize the functions and determine the breaks of the domain:
for i from 1 to n do
   f[i] := args[2*i]; #f[i] is the ith cubic polynomial piece.
   lx[i+1] := args[2*i-1];
end do;
#determine the break points of the domain of the inverse function:
for i from 1 to n do
   ly[i+1] := eval(f[n-i+1], x=lx[n-i+1]);
#print(f[i], lp[i+1], lq[i+1]); #control
end do;
#plot(piecewise(seq( [p<=args[2*i-1], f[i]][], i=1..n) ), p=lp[2]..lp[n+1]);#control
for i from 1 to n do
   L[i] := simplify(evalc([solve(f[i]=y,x)])) assuming y>=ly[n-i+1], y<=ly[n-i+2]: #thanks to Robert Israel for this part
   P[i] := evalf(eval(L[i], y=ly[n-i+1]));
end do;
for i from 1 to n do
   if abs(P[i][1]-lx[i+1])<=0.01 then #suppose that the upper bound of the error is 0.01
      finverse[i] := L[i][1] else
         if abs(P[i][2]-lx[i+1])<=0.01 then
            finverse[i] := L[i][2] else
               finverse[i] := L[i][3]
         end if;   
   end if;
end do;
#for some reason, the following nested loop does not work when i is from 1 to n. It is fine when i is from 1 to n-1, though. I have to do n separately.
#for i from 1 to n-1 do
#   for j from 1 to 3 do
#      if abs(P[i][j]-lp[i+1])<=0.01 then
#         finverse[i] := L[i][j]
#      end if;
#   end do;
#end do;
#concetenate these pieces in a piecewise function.
return piecewise( seq( [y<=ly[i+1], finverse[n-i+1] ][], i=1..n  ) );
end proc:

f := x->piecewise(x <= 2, 13-8*x+4*x^2-x^3, x <= 3, -39.13953960+64.18605280*x-29.08139810*x^2+4.011628300*x^3, x <= 9, 8.013637628-3.472493940*x+.3065739839*x^2-0.1316006154e-1*x^3):
expand(f(x));

f_inverse:=y->Invert(2, 13-8*x+4*x^2-x^3, 3, -39.13953960+64.18605280*x-29.08139810*x^2+4.011628300*x^3, 9, 8.013637628-3.472493940*x+.3065739839*x^2-0.1316006154e-1*x^3);

f_inverse(y);

plot(f(x), x=0..4, y=0..15);

plot(f_inverse(y), y=0..15, discont=true);

g := y->diff(y*f_inverse(y),y):

plot([f_inverse(y), g(y)], y=0..15, discont=true);