Inverse of algebraic expression in Matlab

I need to get inverse expression of the following algebraic expression in matlab.

clc;
clear;
syms r
E(r)= 1-1/2*((1+r)*log2(1+r)+(1-r)*log2(1-r));
g=finverse(E)

Using this command no output is there. I need to see the expression of r in terms of E. Is there any way to get that. Pl somebody help me.

output is

g(r) =
 
Empty sym: 0-by-1

[Thiago Henrique Gomes Lobato]

The problem in your case is that matlab cannot find a close form to your inverse function, so it gives as an output a warning and the empty result, otherwise your code would give the correct result (you could interpret the "r" from g as E and g as "r"). For what reason do you want the inverse expression? If you just want the numeric values I would calculate E in the wished interval and then maybe fit a polynom for the inverse mapping.

[AVM]

@Thiago: Thanks for your reply. Actually I have a function f which has explicitly r dependency but not E. And E is a fucntion of r explicitly. I need to plot that fucntion f with E , but not with r. Then how can I proceed?

Answer 1

[Star Strider]

Possibly:

rfcn = @(E) fsolve(@(r) 1-1/2*((1+r).*log2(1+r)+(1-r).*log2(1-r)) - E, 10);     % Anonymous Function: r(E)
N = 200;
Ev = linspace(0, 1, N);                                                         % Sample  E’ Vector
for k = 1:numel(Ev)
    r(k) = rfcn(Ev(k));
end

figure
plot(Ev, real(r),    Ev, imag(r))
grid

gama = linspace(-1, 1, N);
[Gama,Rm] = ndgrid(gama,r);
G = @(r,gama) atan(r.*tan(N*gama)); %% r is  [0,1];
g = @(r,gama) atan(r.*tan(gama)); %% gama is also [-1,1];
Delta=G(Rm,Gama)-N*g(Rm,Gama);

figure
meshc(gama, Ev, real(Delta))
grid on
xlabel('\gamma')
ylabel('E')
zlabel('\Delta')
view(125,25)

and:

Answer 2

[John D'Errico]

As is often the case, not all simple problems have an algebraic solution. In fact, as you have found, it is trivial to write one where that fails.

However, it is also easy to solve numerically, as long as you have some value that you wish to solve for. If you plot it, you will see the relation is symmetric around the y axis.

fplot(E,[-2,2])

As well, it has singularities at -1 and 1, and does not exceed 1, nor does it have real values for r outside of the open interval (-1,1). That should be clear just by inspection.

pretty(E)
log(1 - r) (r - 1)   log(r + 1) (r + 1)
------------------ - ------------------ + 1
     2 log(2)             2 log(2)

Since you are using a symbolic form, just use vpasolve.

vpasolve(E(r) == 0.5,r,.5)
ans =
0.77994427112328089747637659133002

See that I chose a start point greater than zero, so vpasolve found a positive solution.

I could also have used fzero, of course.

[AVM]

@John: Thanks a lot. Actually I have a function f which has explicitly r dependency but not E. And E is a fucntion of r explicitly. I need to plot that fucntion f with E , but not with r. Then how can I proceed?

[John D'Errico]

A confusing way to write a question, but I fail to see a problem.

You have f, as a function of r. And you have E, also as a function of r.

Now, it seems you wish to plot f, but as a function of E. WTP?

I just showed you how to compute r, given any value of E, as long as E does not exceed 1.

So for ANY value of E <= 1, compute r. Feel free to use a loop, where you now loop over values of E, computing one value of r for every E where you would be interested. Perhaps E will vary from -10 to 1, with some increment. You can save the values of r computed in the loop. That is what god made vectors for.

And given that set of values of r, you can now trivially plot f and E on the same set of axes, since we just said you have the relationship f as a function of r. Again, WTP? You have no choice in how you do this, since we have already established the E( r ) relation lacks a functional inverse.

If you are good at writing code, you can even write a function, call it Einverse. Have it solve the numerical problem using perhaps fzero, for ANY value of E. If you are even better at writing code, you can make Einverse loop over multiple values of E internally, returning one solution r for every value of E provided. Given the numerical function Einverse, now you can trivially plot f(Einverse(E)).

But no matter how you decide to do it, you need to solve the problem of inverting E numerically, and I already showed you exactly how you might do that. So I fail to see wherein lies the problem.

Sigh.

function r = Einverse(E)
  r = nan(size(E));
  E_r = @(r) 1-1/2*((1+r)*log2(1+r)+(1-r)*log2(1-r));
  for Eind = 1:numel(E)
    % note that if E is greater than 1, there is no real solution, so
    % that test will drop around this block of code. As such, r will
    % remain as a NaN
    if E(Eind) == 1
      % special case when E is exactly 1, since we know r there.
      r(Eind) = 0;
    elseif E(Eind) < 1
      r(Eind) = fzero(@(r) E_r(r) - E(Eind),.5);
    end
  end
end

Now, you can trivially plot f(E), since you just do something like this:

E = -10:0.1:1;
plot(E,f(Einverse(E)))

Of course, you need to supply the relationship for f, as a function of r. You never told me that.

[AVM]

@John: Thanks a lot for your detail illustration. Actually I have the follwing function to be plotted.

N=1000; %% Very large number
G=atan(r*tan(N*gama)); %% r is  [0,1];
g=atan(r*tan(gama)); %% gama is also [-1,1];
Dealta=G-N*g;
E= 1-1/2*((1+r)*log2(1+r)+(1-r)*log2(1-r)); %% E is [0,1]

%% E_max=1; when r=0; E_min=0 when r=1

Now I need a mesh plot of the 'Delta ' with gama as x axis and E as Y axis.

I am confused as the Delta has no explicitly E depedency rather it is r dependent.

[John D'Errico]

I'm sorry, but I will not keep chasing a moving target. This question is going far afield in the comments, morphing from one problem to another. I solved your first two questions, but I cannot play the role of long term MATLAB consultant to solve your every problem. If I did, then you will just return with another modified form with an extra variable ot two, but where you are not thinking for yourself.

I'll just suggest that if you think of it along the lines of how I showed you how to solve the first two problems, it will be clear.

[AVM]

Thanks for you reply. Okay, I am trying as you said.