householderVector
householderReflection
roots
This package contains utility functions that are utilized by higher level vector and matrix functions. These functions are usually not useful for an end-user.
Extends from Modelica.Icons.Package (Icon for standard packages).
| Name | Description |
|---|---|
| householderVector
| Calculate a normalized householder vector to reflect vector a onto vector b |
| householderReflection
| Reflect a vector a on a plane with orthogonal vector u |
| roots
| Compute zeros of a polynomial where the highest coefficient is assumed as not to be zero |
Vectors.Utilities.householderVector(a,b);
The function call "householderVector(a, b)" returns the normalized Householder vector
u for Householder reflection of input vector a onto vector b, i.e., Householder vector u is the normal
vector of the reflection plane. Algebraically, the reflection is performed by transformation matrix Q
i.e., vector a is mapped toQ = I - 2*u*u',
with scalar c, |c| = ||a|| / ||b||. Q*a is the reflection of a about the hyperplane orthogonal to u. Q is an orthogonal matrix, i.e.a -> Q*a=c*b
Q = inv(Q) = Q'
a = {2, -4, -2, -1};
b = {1, 0, 0, 0};
u = householderVector(a,b); // {0.837, -0.478, -0.239, -0.119}
// Computation (identity(4) - 2*matrix(u)*transpose(matrix(u)))*a results in
// {-5, 0, 0, 0} = -5*b
| Type | Name | Default | Description |
|---|---|---|---|
| Real | a[:] | Real vector to be reflected | |
| Real | b[size(a, 1)] | Real vector b vector a is mapped onto |
| Type | Name | Description |
|---|---|---|
| Real | u[size(a, 1)] | Housholder vector to map a onto b |
function householderVector
"Calculate a normalized householder vector to reflect vector a onto vector b"
import Modelica.Math.Vectors.norm;
input Real a[:] "Real vector to be reflected";
input Real b[size(a, 1)] "Real vector b vector a is mapped onto";
output Real u[size(a, 1)] "Housholder vector to map a onto b";
protected
Real norm_a=norm(a,2);
Real norm_b=norm(b,2);
Real alpha;
algorithm
assert(norm_b > 0, "Vector b in function housholderVector is zero vector, but at least one element should be different from zero");
assert(norm_a > 0, "Vector a in function housholderVector is zero vector, but at least one element should be different from zero");
alpha := if norm(a + norm_a/norm_b*b,2) > norm(a - norm_a/norm_b*b,2) then norm_a/norm_b else -norm_a/norm_b;
u := (a + alpha*b)/length(a + alpha*b);
end householderVector;
Vectors.Utilities.householderReflection(a,u);
Function "householderReflection(a, u)" performs the reflection of vector
a about a plane orthogonal to vector u (Housholder vector).
Algebraically the operation is defined by
withb=Q*a
where Q is an orthogonal matrix, i.e.Q = I - 2*u*u',
Q = inv(Q) = Q'
a = {2, -4, -2, -1};
u = {0.837, -0.478, -0.239, -0.119};
householderReflection(a,u); // = {-5.0, -0.001, -0.0005, -0.0044}
| Type | Name | Default | Description |
|---|---|---|---|
| Real | a[:] | Real vector a to be reflected | |
| Real | u[size(a, 1)] | householder vector |
| Type | Name | Description |
|---|---|---|
| Real | ra[size(u, 1)] | reflexion of a |
function householderReflection
"Reflect a vector a on a plane with orthogonal vector u"
import Modelica.Math.Vectors;
input Real a[:] "Real vector a to be reflected";
input Real u[size(a, 1)] "householder vector";
output Real ra[size(u, 1)] "reflexion of a";
protected
Real norm_a=Vectors.length(a);
Real h=2*u*a;
algorithm
ra := a - h*u;
// Values close to zero are set to zero.
for i in 1:size(ra, 1) loop
ra[i] := if abs(ra[i]) >= norm_a*1e-12 then ra[i] else 0;
end for;
end householderReflection;
r = Vectors.Utilities.roots(p);
This function computes the roots of a polynomial P of x
P = p[1]*x^n + p[2]*x^(n-1) + ... + p[n-1]*x + p[n+1];
with the coefficient vector p. It is assumed that the first element of p is not zero, i.e., that the polynomial is of order size(p,1)-1.
To compute the roots, the eigenvalues of the corresponding companion matrix C
|-p[2]/p[1] -p[3]/p[1] ... -p[n-2]/p[1] -p[n-1]/p[1] -p[n]/p[1] |
| 1 0 0 0 0 |
| 0 1 ... 0 0 0 |
C = | . . ... . . . |
| . . ... . . . |
| 0 0 ... 0 1 0 |
are calculated. These are the roots of the polynomial.
Since the companion matrix has already Hessenberg form, the transformation to Hessenberg form has not to be performed.
Function eigenvaluesHessenberg
provides efficient eigenvalue computation for those matrices.
r = roots({1,2,3});
// r = [-1.0, 1.41421356237309;
// -1.0, -1.41421356237309]
// which corresponds to the roots: -1.0 +/- j*1.41421356237309
| Type | Name | Default | Description |
|---|---|---|---|
| Real | p[:] | Vector with polynomial coefficients p[1]*x^n + p[2]*x^(n-1) + p[n]*x +p[n-1] |
| Type | Name | Description |
|---|---|---|
| Real | roots[max(0, size(p, 1) - 1), 2] | roots[:,1] and roots[:,2] are the real and imaginary parts of the roots of polynomial p |
encapsulated function roots
"Compute zeros of a polynomial where the highest coefficient is assumed as not to be zero"
import Modelica.Math.Matrices;
input Real p[:]
"Vector with polynomial coefficients p[1]*x^n + p[2]*x^(n-1) + p[n]*x +p[n-1]";
output Real roots[max(0, size(p, 1) - 1),2]=fill(0, max(0, size(p, 1) - 1), 2)
"roots[:,1] and roots[:,2] are the real and imaginary parts of the roots of polynomial p";
protected
Integer np=size(p, 1);
Integer n=size(p, 1) - 1;
Real A[max(n, 0),max(n, 0)] "Companion matrix";
Real ev[max(n, 0),2] "Eigenvalues";
algorithm
if n > 0 then
assert(abs(p[1]) > 0, "Computing the roots of a polynomial with function \"Modelica.Math.Vectors.Utilities.roots\"\n"
+ "failed because the first element of the coefficient vector is zero, but should not be.");
// companion matrix
A[1, :] := -p[2:np]/p[1];
A[2:n, :] := [identity(n - 1),zeros(n - 1)];
// roots are the eigenvalues of the companion matrix
roots := Matrices.Utilities.eigenvaluesHessenberg(A);
end if;
end roots;