using lsqnonlin () for solution of an nonlinear equation

Please consider the following code:

function diff = myfun(D,Y) 

t = 0:30:600; 

l = length(t);

a = 0.5*0.75*10^(-4);

for i = 1:l

Yp(i)=0;

for n = 1:10

y(i,n) = exp(-D*t(i)*(pi^2)*(n^2)/(a^2))/n^2;

Yp(i) = Yp(i) + y(i,n);

end

end

diff = 1- (6*Yp/pi^2)-Y;

%function

Y = [0.0566 0.4432 0.5539 0.6783 0.7303 0.3569 0.4001 0.4278 0.4499 0.4720 ...
    0.4500 0.5157 0.5237 0.5492 0.5590 0.5799 0.5890 0.6000 0.6150 0.6300...
    0.6450];

D0 = 2.0*10^(-12);

options=optimset('LargeScale','on','Display','iter','MaxFunEvals',1e20,'TolFun',2e-50,'TolX',2e-50,'LevenbergMarquardt','on');

[D,resnorm,residual,exitflag,output,lambda]=lsqnonlin(@myfun,D0,[],[],options);

I am trying to get D value using the lsqnonlin() but I am getting the message:

Error using ==> feval
Undefined function or method 'diffusion' for input arguments of type 'double'.

Error in ==> lsqnonlin at 200
            initVals.F = feval(funfcn{3},xCurrent,varargin{:});

Caused by:
    Failure in initial user-supplied objective function evaluation. LSQNONLIN cannot continue.

Please help me with the code why I am getting this error message. What could be the better possible way to obtain the solution for D value.

Answer 1

[Matt J]

I suspect that the code you are running is not the code you have shown. The posted myfun(D,Y) takes 2 arguments, but you are not passing a value for Y to myfun in any way. MATLAB should be complaining that Y is missing, unless your call to lsqnonlin looks something like this,

 [D,resnorm,residual,exitflag,output,lambda]=lsqnonlin(@(D) myfun(D,Y), D0,[],[],options);

Also, if D is a scalar, it would be probably be much simpler to use fminsearch instead of lsqnonlin to minimize norm(diff).

[KB]

Thank you Matt. It is running now but what I have seen the output D is same as the initial D0 value provided for the initialization. So, whatever even if I change the D0 value, the output D value remain same with D0 value with following message 'Optimization terminated: first-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected in trust region model.'This is quite awkward an I have no idea why it is happening.

[Matt J]

It is running now

Now that what? What is your actual code?

[KB]

Sorry to update the actual code:The actual one is as follows

function diff = diffusion(D,Y) 

t = 0:30:600;

l = length(t);

a = 0.5*0.75*10^(-4);

for i = 1:l

Yp(i)=0;

for n = 1:10

y(i,n) = exp(-D*t(i)*(pi^2)*(n^2)/(a^2))/n^2;

Yp(i) = Yp(i) + y(i,n);

end

end

y; 

Yp;

diff = 1- (6*Yp/pi^2)-Y;

%function

Y = [0.05698 0.1700 0.2398 0.2929 0.33510 0.3710 0.4101 0.4256 0.45016 0.47470 0.5010 0.51490 0.5303 0.54908 0.56506 0.5798 0.6012 0.6090 0.6190 0.63650 0.6514]; % experimental data

D0 = 1; 

options=optimset('LargeScale','on','Display','iter','MaxFunEvals',1e20,'TolFun',2e-50,'TolX',2e-50,'LevenbergMarquardt','on');

[D,Resnorm,residual,exitflag,output,lambda]=lsqnonlin(@(D)diffusion(D,Y),D0,[],[],options);

   D,Resnorm

which gives output:

Optimization terminated: first-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected in trust region model.

D =

     1

Resnorm =

       5.8455
even if I change the D0 value, D value return same as D0 with similar Resnorm.
I have passed on the Y values as vector which is experimental data.

[Matt J]

No, this new version of the code still does not run and gives the error below. LevenbergMarquardt is indeed not an optimset option listed in the lsqnonlin documentation,

Error using optimset (line 215)
Unrecognized parameter name 'LevenbergMarquardt'.  Please see the optimset reference page in the
documentation for a list of acceptable option parameters. Link to reference page.

[KB]

I am using 7.12(R2011a) but I don't see that kind of error using the optimset, however the same problem I am facing same problem. Every time, the output D value return same as D0. I am getting following output:

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality   CG-iterations
     0          2         5.84546                             0

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point is less than the selected value of the function tolerance.

criteria details

[Matt J]

Your objective function consists of exponential terms with awkwardly scaled time constants,

      exp(-D*t(i)*(pi^2)*(n^2)/(a^2))

The quantities t(i)*(pi^2)*(n^2)/(a^2) are on the order of 1e+11, so the gradient will be numerically 0 except near D that are scaled very small. When I initialize at D0=0 , I reach a final D=1.4010e-13 after about 90 iterations, and the Resnorm moves from 1.9655 to 2.8453e-04. It looks like 'a' is meant to be much larger, though.

This was on R2013b, running with

 options=optimset('Display','iter','MaxFunEvals',1e20,'TolFun',2e-50,'TolX',2e-50);

[KB]

Matt, you are right. The D value looks nice to me, it has to be within that range because it is a diffusivity value which is an extremely slow process and 'a' value is the diameter of droplet which is in micron range. However, when I am using the same code to solve for D value for another data set of Y, it shows:

                                       Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality   CG-iterations
     0          2         6.66559                     7.02e+008
     1          4         3.50381   8.77753e-009              0            1
Optimization terminated: first-order optimality less than OPTIONS.TolFun,
 and no negative/zero curvature detected in trust region model.

D =

8.7775e-009

Resnorm =

       3.5038
This D value should give similar value in the range of 10^(-12). I tried to do it increasing/decreasing the Tolfun and TOlx, but could not solve the problem.

[Matt J]

If so, why not measure a and other quantities in microns instead of millimeters?

[KB]

The term inside of exponential is unitless, D is in cm^2/s, a in micron or cm, t is s. However, when I am using the same code to solve for D value for another data set of Y.

function diff = diffusion(D,Y)

t = 0:30:600;

l = length(t);

a = 0.5*0.75*10^(-4);

for i = 1:l

Yp(i)=0;

for n = 1:10

y(i,n) = exp(-D*t(i)*(pi^2)*(n^2)/(a^2))/n^2;

Yp(i) = Yp(i) + y(i,n);

end

end

y;

Yp;

diff = 1- (6*Yp/pi^2)-Y;

%function

Y = [0.057854 0.24074 0.33023 0.39484 0.44657 0.49006 0.52771 0.56093 0.59065 0.61749 0.64194 0.66434 0.68496 0.70403 0.72172 0.73818 0.75352 0.76786 0.78127 0.79384 0.80563];

D0 = 0;

options=optimset('LargeScale','on','Display','iter','MaxFunEvals',1e20,'TolFun',2e-50,'TolX',2e-50,'LevenbergMarquardt','on');

[D,Resnorm,residual,exitflag,output,lambda]=lsqnonlin(@(D)diffusion(D,Y),D0,[],[],options);

   D,Resnorm

 it shows:

                                       Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality   CG-iterations
     0          2         6.66559                     7.02e+008
     1          4         3.50381   8.77753e-009              0            1
Optimization terminated: first-order optimality less than OPTIONS.TolFun,
 and no negative/zero curvature detected in trust region model.

D =

8.7775e-009

Resnorm =

       3.5038

This D value supposed to give similar value in the range of 10^(-12). I tried to do it increasing/decreasing the Tolfun and TOlx, but could not solve the problem. The same code worked nice with previous data set but not for this one, although the new Y data set is little different than previous one.

[Matt J]

The finite difference calculations done by lsqnonlin aren't smart enough to know how small D needs to be. You could play with DiffMinChange, I suppose, but scaling the function like below also works. Note also that you can verify by plotting that you've found a valid minimum.

fun=@(D)diffusion(D,Y);

[Dsc,Resnorm]=lsqnonlin(@(Dsc) fun((1e-10)*Dsc),D0,0,[],options);

D=Dsc*1e-10

D_interval=linspace(0,2*D,1000);

plot(D_interval, arrayfun(@(D)norm(fun(D))^2,D_interval));
xlabel 'D'
ylabel 'norm(diff)^2'

[KB]

hello, thanks, I have another problem. Recently, I have installed Matlab 2012a in my personal computer and trying to run the same code. Surprisingly, it is not running with 'LevenbergMarquardt' one, so adopted what you suggested. Still it is not running and shows following error:

Undefined function 'isfinite' for input arguments of type 'diffusion'.

Error in /Applications/MATLAB_R2012a.app/toolbox/shared/optimlib/finDiffEvalAndChkErr.p>finDiffEvalAndChkErr (line 35)

Error in /Applications/MATLAB_R2012a.app/toolbox/shared/optimlib/finitedifferences.p>finitedifferences (line 128)

Error in sfdnls (line 89) [J,~,~,numEvals] = finitedifferences(x,fun,[],lb,ub,valx,[],[], ...

Error in snls (line 178) [A,findiffevals] = sfdnls(xcurr,fvec,Jstr,group,[],funfcn{3},l,u,...

Error in lsqncommon (line 175) [xC,FVAL,LAMBDA,JACOB,EXITFLAG,OUTPUT,msgData]=...

Error in lsqnonlin (line 235) [xCurrent,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...

Error in diffusion_call (line 4) [D,Resnorm,residual,exitflag,output]=lsqnonlin(@(D)diffusion(D,Y),D0,[],[],options);

The same code worked well in past but i have no idea it is not running in my computer.

[Matt J]

Looks like your diffusion function is being confused with one of MATLAB's

http://www.mathworks.com/help/finance/diffusion.html

You might want to rename yours, or at least move it higher on your path.

[KB]

Thanks Matt, I will look into it. Also, sometimes it shows the problem in R2012a: Error using lsqnonlin (line 234) LSQNONLIN requires all values returned by user functions to be of data type double.

Error in diffusion_call (line 4) [D,Resnorm,residual,exitflag,output]=lsqnonlin(@(D)diffusion(D,Y),D0,[],[],options); whereas same code runs well in R2011a.

Scaling and Plotting the values as you suggested gives very different D value than the obtained one using lsqnonlin(). lsqnonlin gives D around 8.9*10^-9 whereas scaling gives 2.8*106(-13). SInce, I know the D value would be in the range of 10^(-13), then scaling works. What if it is totally unknown i.e. at what range the D value should fall.

[Matt J]

LSQNONLIN requires all values returned by user functions to be of data type double.

Make sure your function is returning type double.

lsqnonlin gives D around 8.9*10^-9 whereas scaling gives 2.8*106(-13).

Not for me. lsqnonlin gave me the 2.7*1e(-13) value.

What if it is totally unknown i.e. at what range the D value should fall.

As you can see from my code, I was able to find it by giving nothing nower than a positivity bound D>=0 to lsqnonlin. That's only because the function is unimodal, however. When you have local minima, you would need to be able to provide an accurate initial guess D0.

[KB]

Hello Matt, are you using the same code what I am using in R2011a

function diff = diffusion(D,Y)

t = 0:30:600;

l = length(t);

a = 0.5*0.76*10^(-4);

for i = 1:l

Yp(i)=0;

for n = 1:10

y(i,n) = exp(-D*t(i)*(pi^2)*(n^2)/(a^2))/n^2;

Yp(i) = Yp(i) + y(i,n);

end

end

y;

Yp;

diff = 1- (6*Yp/pi^2)-Y; %% function

Y = [0.0568 0.2376 0.3219 0.4102 0.4487 0.5023 0.5589 0.5787 0.6104 0.6119 0.6418 0.666 0.7015 0.7111 0.7301 0.74109 0.7621 0.7876 0.7901 0.8015 0.81826]; D0 = 0;

options=optimset('Display','iter','MaxFunEvals',1e20,'TolFun',2e-50,'TolX',2e-50);

[D,Resnorm,residual,exitflag,output]=lsqnonlin(@(D)diffusion(D,Y),D0,[],[],options);

D,Resnorm The output is:

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality   CG-iterations
     0          2         6.87754                     7.13e+008
     1          4         3.38556   8.91612e-009              0            0

Local minimum found.

Optimization completed because the size of the gradient is less than the selected value of the function tolerance.

criteria details

D =

8.9161e-009

Resnorm =

    3.3856

[Matt J]

Your code, and its output, are difficult to read. You should use the

button to format it distinctly from your text (like my posted code appears). In any case, the code you've shown is not what I ran. I'm still applying lsqnonlin to the scaled version of the problem:

fun=@(D)diffusion(D,Y);
[Dsc,Resnorm]=lsqnonlin(@(Dsc) fun((1e-10)*Dsc),D0,0,[],options);
 D=Dsc*1e-10;

[KB]

Thanks Matt. You have been really awesome. Could you help me out with this data set:

and I am using similar code:

    function diff = kani(D,Y) 

    t = 0:120:7200;

    l=length(t);
    a = 0.5*0.75*10^(-4); 
    for i = 1:l
        Yp(i)=0;
        for n = 1:10
            y(i,n) = exp(-D*t(i)*(pi^2)*(n^2)/(a^2))/n^2;
            Yp(i) = Yp(i) + y(i,n);
        end
    end
    y; 
    Yp;
    diff = 1- (6*Yp/pi^2)-Y;

Y =    [ 0 0.0035 0.0066 0.0091 0.0121 0.0143 0.0150 0.0161 0.0167 0.0174 0.0188 0.0197 0.0199 0.0210 0.0218 0.0226 0.0229 0.0226 0.0237 0.0247 0.0244 0.0262 0.0249 0.0252 0.0248 0.0267 0.0285 0.0272 0.0279 0.0277 0.0292 0.0294 0.0289 0.0301 0.0282 0.0285 0.0299 0.0285 0.0304 0.0306 0.0309 0.0315 0.0317 0.0310 0.0310 0.0304 0.0339 0.0326 0.0320 0.0323 0.0335 0.0312 0.0332 0.0333 0.0316 0.0311 0.0315 0.0307 0.0313 0.0312 0.0328];

fun=@(D)kani(D,Y);

      D0 = 0;

      options = optimset('Display','iter','MaxFunEvals',1e20,'TolFun',2e-50,'TolX',2e-50);

      [Dsc,Resnorm]=lsqnonlin(@(Dsc) fun((1e-10)*Dsc),D0,0,[],options);

      D=Dsc*1e-10 

      Resnorm

      D_interval=linspace(0,2*D,1000);

      plot(D_interval, arrayfun(@(D)norm(fun(D))^2,D_interval));

      xlabel 'D'

      ylabel 'norm(diff)^2' }

and this is what I got:

I have tried to solve this by scaling the D value very low e^(-20), but still not getting it.