Problem with fsolve when using to find a single variable

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

Answer 1

[Walter Roberson]

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.

[Mike]

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]

Do you have the symbolic toolbox? If so consider using it to create the solution, possibly using a higher Digits setting than the default.

Answer 2

[Matt J]

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

[Matt J]

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]

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.