Linearization of Vector Valued Function

We often come across a vector valued function which needs linearization. Depending on the form though, doing it manually is very tedious and prone to error. In this post, let us try to use a computer algebra system, maxima and try to get a linear approximation of a vector valued function.

Basic Theory

The principle here is that we need to use the Taylor series to obtain linear approximation of the function around a point \bar{x}. Read this post from Khan academy which explains the concept quite well.

Python

We could use sympy available as Python library for symbolic algebra. It is lot more intuitive compared to maxima (atleast for me). Here is the code:

from sympy import *


init_printing( use_unicode=True )
# Define the symbols 
x1, x2, u = symbols( 'x1 x2, u')
a,b,c = symbols( 'a b c')

# Define the function we are going to linearize
f = Matrix( [x1**2 + sin(x2) -1 , -x2**3 + u] )

f_at = f.subs( [(x1,a), (x2,b)] )
J = f.jacobian([x1,x2])
J_at = J.subs( [(x1,a), (x2,b)] )
delta = Matrix([x1,x2]) - Matrix([a,b])

# Tailor Series Linearization
f_approx = (f_at + J_at * delta).expand()

# Getting the coefficient matrix
symbols_list = (x1,x2,u)
coef =  Matrix( [  [rowi_of_f.coeff(s) for s in symbols_list ] for rowi_of_f in f_approx ] )

# A * (x1,x2,u) + b
A = coef
b = f_approx - coef * Matrix([ x1, x2, u] )


Maxima Code

Let us try to get a linear approximation of the example function using maxima. Following is the code:

/* The input function */
f( x1, x2, u ) := [ x1^2 + sin(x2) - 1, -x2^3 + u ];

/* Jacobian of this function */
J(x1, x2, u ) := jacobian( f(x1, x2, u), [x1, x2, u ] );

/* Get jacobian at the point of linearization */
J_at(a,b,c) := at( J(x1,x2,u), [x1=a, x2=b, u=c] );

/* Taylor Series Expansion */
g(x1,x2,u,  a,b,c) := f(a,b,c) + J_at(a,b,c) . ( [x1, x2, u] - [a,b,c] );

/* What we need, */
g( x1,x2,u,   1,0,0 );
(%o8)		[ x2 + 2*(x1-1) ]
		[ u             ]
	

/* can also try, at another linearization point */
g( x1, x2, u,  10, 2, 0 ); 
               [ cos(2) (x2 - 2) + 20 (x1 - 10) + sin(2) + 99 ]
(%o11)         [                                              ]
               [           (- 12 (x2 - 2)) + u - 8            ]


/* g --> A . [x1, x2, u ] + b */
q: list_matrix_entries( g( x1, x2, u, a,b,c ) );
A: coefmatrix( q, [x1, x2, u ] );
b: expand( q - A . [x1,x2,u]  );



Method-2 (Maxima)

The following also gives the same result as above. Following is also less error prone as you don’t need to deal with jacobians.

%(i1) f( x1, x2, u ) := matrix( [ x1^2 + sin(x2) - 1 ], [ -x2^3 + u ] );
%(i2) SS( x1, x2, u ) := taylor( f(x1,x2,u), [x1,x2,u],  [a,b,c], 1 );
                             2
(%o2)/T/ [((- 1) + sin(b) + a ) + (2 a (x1 - a) + cos(b) (x2 - b)) + . . ., 
                                                               3           2
                                                         (c - b ) + ((- 3 b  (x2 - b)) - c + u) + . . .]

q: list_matrix_entries( SS( x1,x2,u ) );
A: coefmatrix( q, [x1,x2,u] );
b: expand( q - A . [x1,x2,u] );

Leave a Reply

Please log in using one of these methods to post your comment:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s