exp*prod*exp(x(lnx)^k) solving equation

I want to solve this equation in MATLAB:

    term2 = for k=1:numel(j) (prod(D(k,:)==exp(P(1).*c.*(log(P(1)).^k/factorial(k))))) end
term1= exp(-(lambda*P(1).*pi.*(r.^2)))

I'm trying to do that:

 function F = root2d(P); 
 lambda = 2*10^-4;
th = -40:-1:-106;
PL1 = 10471285.480509; % (mw)
p1 = 10;
p2 = 6 ; 
p3 = 8 ; 
al = 2.5;
T = 10.^(th./10);
r = (p1*PL1^(-1)./T).^(1/al);
R = (p2*PL1^(-1)./T).^(1/al);
syms P
c = (lambda.*pi.*(R.^2));
j = 1:3;
D = zeros(3,67);
for k = 1:numel(j)
    F(1) = (prod(D(k,:)==exp(P(1).*c.*(log(P(1)).^k/factorial(k))))).*exp(-(lambda*P(1).*pi.*(r.^2)))-P;
end

 %fun = @root2d;   in the command line
%P0 = 0;  in the command line
%P = fsolve(fun,P0)  in the command line

when th = -40

    Error using fsolve (line 269)
FSOLVE requires all values returned by functions to be of data type double
  end.
but when th is a vector (1*67)double (th = -40:-1:-106)
Error in fsolve (line 230)
        fuser = feval(funfcn{3},x,varargin{:});
Caused by: Failure in initial objective function evaluation. FSOLVE cannot continue.
`
Do you have an idea?

[John D'Errico]

This is destined for failure as you are doing it.

1. You seem to be trying to solve a problem for each value of th, which is a vector. But you are starting with a scalar value for P0? Thus, for each value of th, you compute T. Then you compute r and R.

2. As importantly, exponentials are bad things when working in double precision. If you log that expression however, it turns into a sum that will be much better posed, possibly now solvable using double precision arithmetic.

3. You are using syms in conjunction with fsolve. fsolve is not a symbolic solver. Just because you don't know the value of P, does not make this a symbolic problem.

[Marwen Tarhouni]

@John D'Errico:

1 I start with th=40

2 I'm trying to do that:

3 and without syms

endError in fsolve (line 388)
      trustnleqn(funfcn,x,verbosity,gradflag,options,defaultopt,f,JAC,...

Answer 1

[John D'Errico]

Read what I said. Again. TAKE THE LOG!!!!!!! That gives you something like this:

-lambda*pi*r^2*P + lambda*pi*R^2*P*(log(P) + log(P)^2/2 + log(P)^3/6) - log(P) = 0

At least, assuming I did the algebra correctly, always suspect. For a single value of th, better than fsolve is just fzero.

Did you really mean th=40? Or did you intend -40? You originally wrote this:

th = -40:-1:-106;

So 40 is not anywhere in that region. I'll assume -40.

th = -40;

PL1 = 10471285.480509; % (mw)
p1 = 10;
p2 = 6 ; 
p3 = 8 ; 
al = 2.5;
T = 10.^(th./10);
r = (p1*PL1^(-1)./T).^(1/al);
R = (p2*PL1^(-1)./T).^(1/al);

Next, learn what scientific notation does for you, and how to use it. Much easier to write.

lambda = 2e-4;

What are r and R now? ALWAYS check your numbers. do they make sense?

r
r =
       0.1556

R
R =
      0.12684

Next, learn how to write a function. Simplest here is a function handle. Note my careful use of .* and .^ in there.

pfun = @(P) -lambda*pi*r^2*P + lambda*pi*R^2*P.*(log(P) + log(P).^2/2 + log(P).^3/6) - log(P)

Next, before you EVER try to solve a problem, PLOT IT.

ezplot(pfun),grid on,refline(0,0)

To me, it looks like P==1 is a solution, or pretty close to that value.

pfun(1)
ans =
  -1.5212e-05

So P==1 does pretty well.

format long g
[Pval,fval] = fzero(pfun,1)
Pval =
         0.999984788419189
fval =
     -5.23398598767377e-17

0.999984788419189... does much better though.

Now, just loop over the various values for th. Save the results from each pass through your loop in a vector of results.

Now, could I have written the above, without taking the log? Yes. But suppose th was something else, like -106?

r
r =
          67.9203632617184

R
R =
          55.3682121328841

Now, we need to check if there are numerical problems when we use exp. Simpler is to not worry, and work in log form.

[Marwen Tarhouni]

@John D'Errico: thank you for this answer , (th=-40:-1:-106)

you're right if th =-106 P=0 but my problem when i use exp.

i have a problem for th =-40 and many problems for loop iteration th is a vector (67*1)

    th=-40:-1:-106

[John D'Errico]

I fail to see the problem. But then you have not shown what you tried afterwards.

You cannot use the code you have posted. It will fail, and I will not debug a mess of code when I have already shown you how to write this. Using fsolve is silly here, wanting to solve all 67 problems at once. However, I did show how to solve it using fzero, for a unique value of th.

thlist = -40:-1:-106;
Plist = zeros(size(thlist));
for i = 1:numel(thlist)
  th = thlist(i);
  % solve now, for this single instance of th.
  % Save the result in a vector.
  (stuff)
  % thus, once you have computed the value of P...
  Plist(i) = P;
end

I suppose, you could use arrayfun, if you understand how to use that tool. But simply not worth the effort to explain it, when a trivial loop will be as fast and good.

So WTP?

[Marwen Tarhouni]

@John D'Errico: Thank you for accuracy & for your detailed answer