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

1 次查看(过去 30 天)
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 个评论
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 CenterFile Exchange 中查找有关 Calculus 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by