In this post, I share a numerical Jacobian matrix calculation method with matlab code.

Actually there is a function in Matlab inherently, but it is very complex. ( look at the function, NumJac ), So I made a very simple version. I wish this can help you.

You can download from here.

The attachment file includes

1. NumJacob.m : main file it generates jacobian matrix.

2. Demo.m : demo file

3. test_func.m : test function to show its demo

The below is a simple example.

>> NumJacob(@cos,1)

ans =

-0.8415

>> x0=[0;1]; other_param=[1;2];

>> df=NumJacob(@test_func,x0,other_param)

df =

1.0000 0

0 2.0000

Because the code is very easy, probably you can easily understand the codes.

Then, good luck ^^

if you have any question, please leave me a comment below.

Added on Feb 07. 2014

I am uploading another demo file to give a solution of Mahmudul’s request (for detail, see the comment).

This demo file is to get a numerical Jacobian for the below function

y1′ = y2*y3; y1(0) = 0

y2′ = -y1*y3; y2(0) = 1

y3′ = -0.51*y1*y3; y3(0) =1

Download this demo file.

—————————————————————————————————

I am Youngmok Yun, and writing about robotics theories and my research.

My main site is http://youngmok.com, and Korean ver. is http://yunyoungmok.tistory.com.

—————————————————————————————————

### Like this:

Like Loading...

MahmudulDear Youngmok, I am a post graduate research student at University and using Matlab for my modelling purpose. I just wonder if you could clarify what the 2nd and 3rd input arguments of the “function df=NumJacob(f,x0,varargin)”. Actually, I would be grateful to you, if you just let me know the suitable command for evaluating the jacobian matrix of the following system od ODE:

y1′ = y2*y3; y1(0) = 0

y2′ = -y1*y3; y2(0) = 1

y3′ = -0.51*y1*y3; y3(0) =1

Many thanks in advance!

adminPost authorI am uploading another example to solve your problem. Read the post one more carefully, and then you can download the demo file for you.

MahmudulDear Youngmok,

Your help is much appreciated. I found your uploaded files quite handy.

Thanks again!

adminPost authorIt is my plesure. I am happy to help someone who needs my help. Good luck to your research.

MahmudulHi Youngmok,

Your previous file was helpful to me. However, I am still struggling to find the way of my problem.

My problem is something like given below:

The system is in the form of y’ = f(time,y) and is described below:

y’ = f(time,y)

if time > 10

y1′ = y2*y3; y1(0) = 0

y2′ = -y1*y3; y2(0) = 1

y3′ = -0.51*y1*y3; y3(0) =1

else

y1′ = y2; y1(0) = 0

y2′ = -y1*y3; y2(0) = 1

y3′ = -0.51*y1*y3; y3(0) =1

end

Now, actually I want to evaluate the jacobian matrix at some fixed points ( i.e. where y’ = 0) at different time stage…

So, it needed to find the fixed points ( i.e. values of all variables by setting y’=0) and then evaluate the jacobian matrix at these points rather than at initial conditions, as your previously uploaded file ( where 0 1 1 were the initial values of the system). Please note that, my real system of odes is nonlinear ( i.e involves ^2 and ^4 terms). I’ll be extremely grateful to you, if you could give me a way out.

Millions of Thanks in advance…

adminPost authorDear Mahmudul,

Sorry for late reply.

In this kind of case, I made the function to receive other parameters.

For example,

>> x0=[0;1]; other_param=[1;2];

>> df=NumJacob(@test_func,x0,other_param)

df =

1.0000 0

0 2.0000

Here, “other_param” is the other parameters.

It means that you can send your another variable (e.g., time) using this method.

Then, Good luck

-Mok-

MahmudulDear Youngmok,

Many thanks for your reply. Are you happy to let me know you skype ID, if you have? And possibly with the information about the best time to ring you…..please

burashkaHello Youngmok.

Would the code work with any system? my task is given an .m file with equations, the function can give the numerical jacobian given the initial guesses (numerical jacobian would be used to solve a system of equations using newton raphson) Would your code work with lets say 32 equations? do i need to modify anything?

Thank you

adminPost authorThere is no limitation for the dimension. Yes. you can use the Numerical jacobian function without any change.

BurashkaThank you for the reply. my test file is given in the following format:

function f = test(x)

f = [

f(1)=Pg1-Pl(1)-real((V1*cos(d1)+1.*i*V1*sin(d1))*((-2.8289-28.289*i)*(V2*cos(d2)-1.*i*V2*sin(d2))+(2.8289+28.289*i)*(V1*cos(d1)-1.*i*V1*sin(d1))));

f(2)=Qg1-Ql(1)-imag((V1*cos(d1)+1.*i*V1*sin(d1))*((-2.8289-28.289*i)*(V2*cos(d2)-1.*i*V2*sin(d2))+(2.8289+28.289*i)*(V1*cos(d1)-1.*i*V1*sin(d1))))];

etc.. i have 32 equations in this format.

and my initial guess is

x0 =[ 1 0 1 0 0 0 1 0 1 0 1 0 1 0 1 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0];

when i run your jacobian on my test file, it is giving me an error. is it because my unknowns are not in terms of x?

youngmok yunx should be a column vector, try again.

Burashkaas in

f = [ Pg1-Pl(1)-real((V1*cos(d1)+1.*i*V1*sin(d1))*((-2.8289-28.289*i)*(V2*cos(d2)-1.*i*V2*sin(d2))+(2.8289+28.289*i)*(V1*cos(d1)-1.*i*V1*sin(d1))));

Qg1-Ql(1)-imag((V1*cos(d1)+1.*i*V1*sin(d1))*((-2.8289-28.289*i)*(V2*cos(d2)-1.*i*V2*sin(d2))+(2.8289+28.289*i)*(V1*cos(d1)-1.*i*V1*sin(d1))))];

it still doesnt work. my unknowns are v1, v2, d1, d2…etc do they need to be in terms of x like x(1), x(2)…

Burashkai just understood what you meant. let me try x as a column vector

youngmok yunif you want to get a partial derivative with respect to v1,v2,d1,d2, you have to include those in x. If those variables (v1,v2,d1,d2) are just parameters of the function, you don’t need to add it in x. If v1,v2,d1,d2 are just parameters but variables, those can be transmitted as ‘other_param’ in NumJacob(@test_func,x0,other_param),

Good luck.

Burashkaim sorry i dont understand. lets say v1 and d1 are just variables (v1 = 1, d1 = 0) then other param would be = [1;0]?

and if i need the partial derivatives w.r.t to v2, d2 then x just becomes [0;0]? (if v2 = 0 and d2= 0)

thanks for all the help. much appriciated