Jacobian of a Matrix

[salvomic]

hi all,
here there is a CAS program to calc the Jacobian of a Matrix.

Input: [f(1), f(2), ...], [x,y,...]

Enjoy!

Salvo Micciché

Code:


#cas
jacob(args):=
// Jacobian Matrix by Salvo Micciché
// input vectorial expression, vector of variables
BEGIN
local argv, argc, mat, f, var, fn, j, k,  gr, vd;
argv:=[args];
argc:=size(argv);
IF argc !=2 THEN
return "Input:[f1(x),f1(y),f1(z)...], [x,y,z,...]"; 
ELSE
f:=argv(1);
var:=argv(2);
fn:=size(f);
vd:=size(var);
mat:=makemat(0,fn,vd);
FOR j FROM 1 TO fn DO // gradients
gr:=grad(f(j),var);
FOR k FROM 1 TO vd DO // items
mat[j,k]:=gr(k);
END; // for k
END; // for j
return mat;
END; // if-else

END;
#end

[Rudi]

hi all,
I've improved the Code of the Jacobian Matrix by salvomic - 15.05.2015 21:20

Enjoy!

Rudi Steeger

For Example:

jacob(grad([e^(x*y^2)*sin(z)],[x,y,z]),[x,y,z]);

The same Example with the new Feature in this Code:
jacob2([e^(x*y^2)*sin(z)],[x,y,z])
==>
[[[y^4*e^(x*y^2)*sin(z)],[2*y*(x*y^2+1)*e^(x*y^2)*sin(z)],[y^2*cos(z)*e^(x*y^2)]],[[2*y*(x*y^2+1)*e^(x*y^2)*sin(z)],[2*x*(2*x*y^2+1)*e^(x*y^2)*sin(z)],[2*x*y*cos(z)*e^(x*y^2)]],[[y^2*cos(z)*e^(x*y^2)],[2*x*y*cos(z)*e^(x*y^2)],[-e^(x*y^2)*sin(z)]]];

Code:
#cas
jacob2(args):=
// Jacobian Matrix by Salvo Micciché
// input vectorial expression, vector of variables
BEGIN
local argv, argc, mat, f, var, fn, fg, j, k, gr, vd;
argv:=args;
argc:=size(argv);
IF argc !=2 THEN
return "Input:[f1(x),f1(y),f1(z)...], [x,y,z,...]";
ELSE
f:=argv(1);
var:=argv(2);
fn:=size(f);
vd:=size(var);
IF fn:=1 THEN
fg:=grad(f(1),var);
f:=fg;
fn:=size(f);
END;
mat:=makemat(0,fn,vd);
FOR j FROM 1 TO fn DO // gradients
gr:=grad(f(j),var);
FOR k FROM 1 TO vd DO // items
mat[j,k]:=factor(gr(k));
END; // for k
END; // for j
return mat;
END; // if-else
END;
#end

[Han]

A slightly shorter program:

Code:


#cas
jacob(args):=
begin
local argv, argc, mat, f, var, fn, j, k,  gr, vd;
argv:=[args];
argc:=size(argv);
IF argc !=2 THEN
return "Input:[f1(x1,...,xn),...,fm(x1,...,xn)], [x1,...,xn]"; 
ELSE
return transpose(diff(argv[1],argv[2]));
END;
end;
#end

Basically, the Jacobian is: transpose(diff([f1,f2,...,fm], [x1,x2,...,xn]))

[Rudi]

hi all,
Wow very nice, the Jacobi matrix contains only the first derivatives. However, the jaboc function calculates the 2nd derivatives. Corresponds essentially to the Hessian matrix. For me it was important to understand the Jacob function in connection with matrices. I am therefore able to write similar functions.
Thank you very much.

P.S.
Translated with Google.

Enjoy!

Rudi Steeger

[sitomix]

Can you make a new function to evaluate jacobian in one point, for example: jacob( [x*y,y*y],[x,y],[x=2,y=5]) ?

[salvomic]

(04-09-2018 12:49 PM)sitomix Wrote:  Can you make a new function to evaluate jacobian in one point, for example: jacob( [x*y,y*y],[x,y],[x=2,y=5]) ?

I'll thing about it, thanks. As soon as I'll have some spare time...

Salvo

[Arno K]

I improved Han's program a little bit, now you can enter what you desire, with and without substitution:

Code:

#cas
 jacob(args):=
 begin
 local argv,argc,mat,f;
 LOCAL var,fn,j,k,gr,vd;
 argv:=[args];
 argc:=size(argv);
 IF argc == 3 THEN
   return subst(transpose(diff(argv[1],argv[2])),argv[3]);
  END;
 IF argc == 2 THEN
   return transpose(diff(argv[1],argv[2]));
  END;
 return "Input:[f1(x1,...,xn),...,fm(x1,...,xn)], [x1,...,xn][,[x1=a1,...xn=an]]"; 
 end;
#end

For its usage see the following picture.
Hope that helps
Arno

[Arno K]

(04-09-2018 12:49 PM)sitomix Wrote:  Can you make a new function to evaluate jacobian in one point, for example: jacob( [x*y,y*y],[x,y],[x=2,y=5]) ?

I am deeply impressed how people wanting our ( that is everybody really involved in things like improving programs, for example) help, finally appreciate this, a tiny "thank you" by sitomix would have been nice, this is one thing my parents taught me, when I was young.
Arno