RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
newton_raphson.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
3  * 2020 Freie Universität Berlin
4  *
5  * This file is subject to the terms and conditions of the GNU Lesser General
6  * Public License v2.1. See the file LICENSE in the top level directory for more
7  * details.
8  */
9 
23 #include "vector.h"
24 #include "matrix.h"
26 
27 //n is the size of x0 that is equal to the column number of J
28 uint8_t newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[],
29  double eps, uint8_t max_it_num, vector_t est_x_arr[],
30  void (*get_non_lin_sys)(vector_t x_arr[],
31  vector_t f_vec[]),
32  void (*get_jacobian)(vector_t x_arr[], matrix_t J[][n]))
33 {
34  uint8_t iter_num;
35  double step;
36  vector_t prev_x_arr[n];
37  matrix_t J[f_length][n];
38  matrix_t pinv_J[n][f_length];
39  vector_t delta_x[n];
40  vector_t f_vec[f_length];
41 
42  vector_copy(n, x0_arr, prev_x_arr);
43  step = eps;
44  iter_num = 0;
45 
46  while ((step >= eps) && (iter_num < max_it_num)) {
47  // compute the Jacobian matrix for the prev_x_arr.
48  get_jacobian(prev_x_arr, J);
49  // compute the inverse of J (J^-1)
50  moore_penrose_get_pinv(f_length, n, J, pinv_J);
51  // compute the value of f(x)
52  get_non_lin_sys(prev_x_arr, f_vec);
53  //delta_x = J1^-1*f(x)
54  matrix_mul_vec(n, f_length, pinv_J, f_vec, delta_x);
55  //improve x: x = x - delta_x
56  vector_sub(n, prev_x_arr, delta_x, est_x_arr);
57  //next step
58  step = vector_get_euclidean_distance(n, est_x_arr, prev_x_arr);
59  //update prev_x
60  vector_copy(n, est_x_arr, prev_x_arr);
61  iter_num++;
62  }
63 
64  return iter_num;
65 }
moore_penrose_get_pinv
int8_t moore_penrose_get_pinv(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t pinv_A[n][m])
Calculate the Moore–Penrose inverse of a rectangular matrix.
Definition: moore_penrose_pseudo_inverse.c:38
matrix_mul_vec
void matrix_mul_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n], matrix_t vec[n], matrix_t dst_arr[m])
Compute the multiplication of a matrix with a column vector.
Definition: matrix.c:434
vector_get_euclidean_distance
vector_t vector_get_euclidean_distance(uint8_t length, vector_t vec1[], vector_t vec2[])
Compute the Euclidean distance between two vectors.
Definition: vector.c:163
vector_t
#define vector_t
Define the data type of the vector elements.
Definition: vector.h:33
matrix.h
Matrix computations.
moore_penrose_pseudo_inverse.h
Moore–Penrose algorithm to compute the pseudo-inverse of a matrix.
vector_copy
void vector_copy(uint8_t size, vector_t src_arr[], vector_t dest_arr[])
Copy the elements of the source vector to the destination vector.
Definition: vector.c:37
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
vector.h
Vector computations.
newton_raphson
uint8_t newton_raphson(uint8_t f_length, uint8_t n, vector_t x0_arr[], double eps, uint8_t max_it_num, vector_t est_x_arr[], void(*get_non_lin_sys)(vector_t x_arr[], vector_t f_vec[]), void(*get_jacobian)(vector_t x_arr[], matrix_t J[][n]))
Implements the Newton–Raphson algorithm.
Definition: newton_raphson.c:28
vector_sub
void vector_sub(uint8_t size, vector_t a_vec[], vector_t b_vec[], vector_t a_minus_b[])
Compute the subtraction of two vectors.
Definition: vector.c:94