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);