Prev: Close figure before show it
Next: how to read an Xls file that is protected with password
From: derek on 7 Jun 2010 11:00
(version R2009b)

I've written program where I am solving a PDE such that for every timestep I am solving a system of nonlinear equations F(X) = 0 with 100+ degrees of freedom. I can use fsolve for it, but I find that all the algorithms built into fsolve have trouble converging for one reason or another. However, if I write a simple routine that does just basic Newton's method I find it works much better. This would be fine, but to use Newton's method I either need to calculate the actual Jacobian (I'd rather avoid that if possible as it would be very tedious), or I can use a finite difference approximation of the Jacobain. This is very slow though, since if I have N degrees of freedom, I need N evaluations of the function F(X) in order to estimate the Jacobian.

In the fsolve function, when you use the trust-region-reflective algorithm, you can pass in a JacobPattern sparse matrix that has the pattern for the Jacobian so that it can more quickly compute the finite difference approximation of the Jacobian. Is there a way I could extract this part of fsolve so that I can more quickly get the finite difference Jacobain estimate for my own solver?

Looking through the code for fsolve, starting on line 355 for the trust-region-reflective algorithm it has the following code:

% Execute algorithm
if algorithmflag == 1 % trust-region reflective
if ~gradflag
Jstr = optimget(options,'JacobPattern',defaultopt,'fast');
if ischar(Jstr)
% options.JacobPattern is the default: 'sparse(ones(jrows,jcols))'
Jstr = sparse(ones(Jrows,Jcols));
end
checkoptionsize('JacobPattern', size(Jstr), Jcols, Jrows);
else
Jstr = [];
end
computeLambda = 0;
% Set MaxFunEvals appropriately for trust-region-reflective
defaultopt.MaxFunEvals = '100*numberOfVariables';

[x,FVAL,LAMBDA,JACOB,EXITFLAG,OUTPUT,msgData]=...
snls(funfcn,x,LB,UB,verbosity,options,defaultopt,f,JAC,caller,...
Jstr,computeLambda,mtxmpy,detailedExitMsg,optionFeedback,varargin{:});

So it uses this snls function, located at *\MATLAB\2009b\toolbox\optim\optim\private\snls.m. This I think is the actual solver itself, since looking at the code for it it stands for 'sparse nonlinear least squares' solver. In that function on lines 183 and 201 it uses the function finitedifferences to evaluate the Jacobian.

This finitedifferences function, located at *\MATLAB\2009b\toolbox\shared\optimlib\finitedifferenes.m, seems to be the function I would want to use, but I can't figure out how to use it to quickly evaluate the Jacobian from the JacobPattern sparse matrix, instead of one row at a time like I am already doing. Would someone know how I could do that?
From: Paul Kerr-Delworth on 8 Jun 2010 06:27
"derek " <derek(a)che.utexas.edu> wrote in message <huj1i6$q0e$1(a)fred.mathworks.com>...
> (version R2009b)
>
> I've written program where I am solving a PDE such that for every timestep I am solving a system of nonlinear equations F(X) = 0 with 100+ degrees of freedom. I can use fsolve for it, but I find that all the algorithms built into fsolve have trouble converging for one reason or another. However, if I write a simple routine that does just basic Newton's method I find it works much better. This would be fine, but to use Newton's method I either need to calculate the actual Jacobian (I'd rather avoid that if possible as it would be very tedious), or I can use a finite difference approximation of the Jacobain. This is very slow though, since if I have N degrees of freedom, I need N evaluations of the function F(X) in order to estimate the Jacobian.
>
> In the fsolve function, when you use the trust-region-reflective algorithm, you can pass in a JacobPattern sparse matrix that has the pattern for the Jacobian so that it can more quickly compute the finite difference approximation of the Jacobian. Is there a way I could extract this part of fsolve so that I can more quickly get the finite difference Jacobain estimate for my own solver?
>
> Looking through the code for fsolve, starting on line 355 for the trust-region-reflective algorithm it has the following code:
>
> % Execute algorithm
> if algorithmflag == 1 % trust-region reflective
> if ~gradflag
> Jstr = optimget(options,'JacobPattern',defaultopt,'fast');
> if ischar(Jstr)
> % options.JacobPattern is the default: 'sparse(ones(jrows,jcols))'
> Jstr = sparse(ones(Jrows,Jcols));
> end
> checkoptionsize('JacobPattern', size(Jstr), Jcols, Jrows);
> else
> Jstr = [];
> end
> computeLambda = 0;
> % Set MaxFunEvals appropriately for trust-region-reflective
> defaultopt.MaxFunEvals = '100*numberOfVariables';
>
> [x,FVAL,LAMBDA,JACOB,EXITFLAG,OUTPUT,msgData]=...
> snls(funfcn,x,LB,UB,verbosity,options,defaultopt,f,JAC,caller,...
> Jstr,computeLambda,mtxmpy,detailedExitMsg,optionFeedback,varargin{:});
>
> So it uses this snls function, located at *\MATLAB\2009b\toolbox\optim\optim\private\snls.m. This I think is the actual solver itself, since looking at the code for it it stands for 'sparse nonlinear least squares' solver. In that function on lines 183 and 201 it uses the function finitedifferences to evaluate the Jacobian.
>
> This finitedifferences function, located at *\MATLAB\2009b\toolbox\shared
> \optimlib\finitedifferenes.m, seems to be the function I would want to use, but I > can't figure out how to use it to quickly evaluate the Jacobian from the
> JacobPattern sparse matrix, instead of one row at a time like I am already doing.
> Would someone know how I could do that?

Hi Derek,

As you're solving a PDE, I wonder if you could use the Partial Differential Equation toolbox instead of writing your own solver. The documentation for the toolbox is available online

<<http://www.mathworks.com/access/helpdesk/help/toolbox/pde/>>

Best regards,

Paul
From: derek on 8 Jun 2010 09:46
> Hi Derek,
>
> As you're solving a PDE, I wonder if you could use the Partial Differential Equation toolbox instead of writing your own solver. The documentation for the toolbox is available online
>
> <<http://www.mathworks.com/access/helpdesk/help/toolbox/pde/>>
>
> Best regards,
>
> Paul

The PDE toolbox is basically just a simplified and stripped-down COMSOL, which isn't what I need for this problem (besides I have COMSOL anyway). My entire problem is 8 coupled equations, some of which are nonlinear PDE's, some are nonlinear ODE's, and others are just equations that relate certain variables. Even though its 1D, the domain includes a moving mesh and changing topology (i.e. I have to merge elements partway through the simulation) so it's really just too complicated for something like the PDE toolbox or even COMSOL. I think my only recourse then is to analytically work out the Jacobian if I want my solver to go any faster.
From: Paul Kerr-Delworth on 11 Jun 2010 08:00
"derek " <derek(a)che.utexas.edu> wrote in message <hulhis$6m2$1(a)fred.mathworks.com>...
> > Hi Derek,
> >
> > As you're solving a PDE, I wonder if you could use the Partial Differential Equation toolbox instead of writing your own solver. The documentation for the toolbox is available online
> >
> > <<http://www.mathworks.com/access/helpdesk/help/toolbox/pde/>>
> >
> > Best regards,
> >
> > Paul
>
> The PDE toolbox is basically just a simplified and stripped-down COMSOL, which isn't what I need for this problem (besides I have COMSOL anyway). My entire problem is 8 coupled equations, some of which are nonlinear PDE's, some are nonlinear ODE's, and others are just equations that relate certain variables. Even though its 1D, the domain includes a moving mesh and changing topology (i.e. I have to merge elements partway through the simulation) so it's really just too complicated for something like the PDE toolbox or even COMSOL. I think my only recourse then is to analytically work out the Jacobian if I want my solver to go any faster.

Hi Derek,

There is a function available that computes a numerical computes a Jacobian. The function is called numjac - type help numjac at the MATLAB command prompt for more information.

If your Jacobian is sparse, numjac allows you to specify a sparsity pattern.

Hope this helps

Best regards,

Paul
 | 
Pages: 1
Prev: Close figure before show it
Next: how to read an Xls file that is protected with password