Does anybody know how to obtain an adecuate step for finite differences in Matlab?

Hi, I have a vector of dim Nx1,fnl0 which depends on a vector X0 of dim 3Nx1. I would like to obtain the parcial derivatives of this vectorfnl0 respective X0 by finite differences, in order to calculate the Jacobian, dfnl0_dX . I´m trying to construct the Jacobian column by column, by concatenating the results of each step k:

X0(vector dim 3Nx1) fnl0(vector dim Nx1) delta_X0=sqrt(eps)*norm(X0) dfnl0_dX=Zeros[N,3N] % Jacobian

for k=1:3*N

X0pert=X0; X0pert(k,1)=X0pert(k,1)+delta_X0; fnl0_pert=f_nonlinear(dx0pert,param); derivative= (fnl0_pert - fnl0) / delta_X0; dfnl0_dX=horzcat(dfnl0_dX,derivative);

end

The Problem is that the elements of vector X0 are very different between them, some elements 1e-20 and others of order 1e-6, and so it seems that I am not obtaining correct results, as I have been trying to do the same Operation with different delta_X0 and the results vary a lot.

Could anyone help me in obatining an adecuate step delta_X0 to solve this problem?

Thank your very much!

Answer 1

[Torsten]

I looked into a sophisticated solver.

There,

Delta_X0(k) = sqrt(Uround*max(1e-5,abs(X0(k))))

where Uround is the SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0.

Note that Delta_X0 is a vector with different values depending on the component.

Best wishes

Torsten.

[student]

thank you, Torsten.I have a question, what does 1.0D0 mean? and is 1e-5 a fixed value?

[Torsten]

1.0D0 means "1 in double precision" and yes, 1e-5 is a fixed value.

Best wishes

Torsten.

[student]

Oh, thank you. What I don#t see know is how to define Uround..

[Torsten]

The default value in the code is 1e-16.

Best wishes

Torsten.

[student]

Hi Torsten, I'm wondering if you could tell me where exactly did you found the formula of the step for finite differences

Delta_X0(k) = sqrt(Uround*max(1e-5,abs(X0(k)))) where Uround is the SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0.

Thank you very much

[Torsten]

http://www.unige.ch/~hairer/prog/stiff/radau5.f

Search for the line

DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE)))

Best wishes

Torsten.