Problem with fsolve when using to find a single variable

1 次查看（过去 30 天）
Mike
Mike 2013-2-18
Hello,
I need help with a problem I am trying to solve using fsolve. The code is as follows:
C=0.25; %Nitrogen Concentration after time t in at%
Co=0.5; %Initial Nitrogen Concentration in at%
D=1.6613e-08; %Diffusion Coefficient at 450C in cm^2/sec
x=0.003; %Penetration depth in cm
t0=1; %Initial guess for time in sec
fun=@(t) [Co*erfc(x/(2*(D*t)^0.5))-C];
time=fsolve(fun,t0)
I would think that this would be a fairly straight forward application of fsolve but I am unable to solve using Matlab. The answer is time equals about 600 seconds which I found using Excel.
Can someone please let me know what I am doing wrong?
Thank you very much,
Mike

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

回答（2 个）

Walter Roberson
Walter Roberson 2013-2-18
The calculation is sensitive to round-off. If not enough guard digits are used, the solution of roughly 595.40673113197968463 is only one of at least 14 possible solutions, the other 13 of which are imaginary.
  2 个评论
Mike
Mike 2013-2-18
Thank you for the suggestion...how do I incorporate guarde digits in the code? Also can I tell Matlab to ignore imaginary answers?
Mike
Walter Roberson
Walter Roberson 2013-2-18
Do you have the symbolic toolbox? If so consider using it to create the solution, possibly using a higher Digits setting than the default.

请先登录，再进行评论。

Matt J
Matt J 2013-2-18
编辑：Matt J 2013-2-18
I find I get pretty good results if I make the change of variables
z=x/(2*(D*t)^0.5);
It also avoids possibly illegal non-differentiabilities from the square roots you're taking
fun=@(z) Co*erfc(z)-C;
z=fzero(fun,t0);
t=(x/2/z/sqrt(D))^2
  2 个评论
Matt J
Matt J 2013-2-18
In fact, you don't even have to solve iteraitvely. Just use ERFCINV,
z=erfcinv(C/Co);
t=(x/2/z/sqrt(D))^2
Walter Roberson
Walter Roberson 2013-2-18
Technically, Mike does not have a sqrt() call: make has a ^0.5 call, which MATLAB defines in terms of logs. On the other hand, sqrt() and ^0.5 both have problems with differentiation at 0.

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by