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.

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.

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

 

 

 

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

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

      Reply
  2. 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…

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

      Reply
  3. 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

    Reply
  4. 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

    Reply
    1. admin Post author

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

      Reply
  5. 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?

    Reply
  6. 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)…

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

    Reply
  8. 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

    Reply
  9. Pingback: Fix Comsol Error Failed To Evaluate Variable Jacobian Windows XP, Vista, 7, 8 [Solved]

Leave a Reply