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

1 次查看（过去 30 天）
Marwen Tarhouni
Marwen Tarhouni 2018-6-16
评论： Marwen Tarhouni 2018-6-18
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?
  2 个评论
John D'Errico
John D'Errico 2018-6-16
编辑：John D'Errico 2018-6-16
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
Marwen Tarhouni 2018-6-16
@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,...

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2018-6-16
编辑：John D'Errico 2018-6-16
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.
  3 个评论
显示 1更早的评论
John D'Errico
John D'Errico 2018-6-17
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
Marwen Tarhouni 2018-6-18
@John D'Errico: Thank you for accuracy & for your detailed answer

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by