# Numerical Jacobian matrix calculation method with matlab code

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.

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.

Added on Mar/21 2015

If you are interested in calculating numerical Jocobian matrix for a complex number, please read this post.

http://youngmok.com/matlab-code-for-numerical-calculation-of-jacobian-matrix-for-a-complex-number/

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

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.

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

## 19 thoughts on “Numerical Jacobian matrix calculation method with matlab code”

1. Mahmudul

Dear 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!

1. admin Post author

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

2. Mahmudul

Dear Youngmok,

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

Thanks again!

1. admin Post author

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

3. Mahmudul

Hi 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…

1. admin Post author

Dear 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-

4. Mahmudul

Dear 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

5. burashka

Hello 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

1. admin Post author

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

6. Burashka

Thank 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?

7. Burashka

as 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)…

8. youngmok yun

if 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.

9. Burashka

im 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