RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
fsolve_test.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 <stdio.h>
24 #include <stdlib.h>
25 #include <math.h>
26 #include <complex.h>
27 #include <float.h>
28 
29 #include "vector.h"
30 #include "matrix.h"
31 #include "damped_newton_raphson.h"
32 #include "newton_raphson.h"
33 #include "fsolve.h"
34 
35 /********************** f1 ******************************/
36 
37 /*
38  * f1(x1,x2) = [x1^3 + x2 - 1;
39  * x2^3 - x1 + 1]
40  */
54 void get_non_lin_sys_f1_(vector_t x_arr[], vector_t f1_vec[])
55 {
56 
57  f1_vec[0] = pow(x_arr[0], 3) + x_arr[1] - 1;
58  f1_vec[1] = pow(x_arr[1], 3) - x_arr[0] + 1;
59 
60 }
61 
62 /*
63  * J1(x1,x2) = [3*x1^2, 1;
64  * -1, 3*x2^2]
65  */
66 
80 void get_non_lin_sys_J1_(vector_t x_arr[], matrix_t J1[][2])
81 {
82  J1[0][0] = 3 * pow(x_arr[0], 2);
83  J1[0][1] = 1;
84  J1[1][0] = -1;
85  J1[1][1] = 3 * pow(x_arr[1], 2);
86 }
87 
88 /****************************** f2 ********************************/
89 
90 /*
91  * f2(x1,x2,x3) = [3*x1-cos(x2*x3)-1/2;
92  * x1^2-81*(x2+0.1)^2+sin(x3)+1.06;
93  * exp(-x1*x2)+20*x3+(10*pi-3)/3
94  * ]
95  */
110 void get_non_lin_sys_f2_(vector_t x_arr[], vector_t f2_vec[])
111 {
112 
113  f2_vec[0] = 3 * x_arr[0] - cos(x_arr[1] * x_arr[2]) - 0.5;
114  f2_vec[1] = pow(x_arr[0], 2) - 81 * pow((x_arr[1] + 0.1), 2)
115  + sin(x_arr[2]) + 1.06;
116  f2_vec[2] = exp(-x_arr[0] * x_arr[1]) + 20 * x_arr[2]
117  + (10 * M_PI - 3) / 3;
118 }
119 
120 /*
121  * [ 3, x3*sin(x2*x3), x2*sin(x2*x3)]
122  * Jf2 = [ 2*x1, -162*x2 - 81/5, cos(x3)]
123  * [ -x2*exp(-x1*x2), -x1*exp(-x1*x2), 20]
124  */
139 void get_non_lin_sys_J2_(vector_t x_arr[], matrix_t J2[][3])
140 {
141  J2[0][0] = 3;
142  J2[0][1] = x_arr[2] * sin(x_arr[1] * x_arr[2]);
143  J2[0][2] = x_arr[1] * sin(x_arr[1] * x_arr[2]);
144  J2[1][0] = 2 * x_arr[0];
145  J2[1][1] = -162 * x_arr[1] - 16.2;
146  J2[1][2] = cos(x_arr[2]);
147  J2[2][0] = -x_arr[1] * exp(-x_arr[0] * x_arr[1]);
148  J2[2][1] = -x_arr[0] * exp(-x_arr[0] * x_arr[1]);
149  J2[2][2] = 20;
150 }
151 
152 /****************************** f3 ********************************/
153 /*
154  *
155  * f3(x1,x2,x3) = [exp(-x1*x2) + log(x1)-exp(-2);
156  * exp(x1)-sqrt(x3)/x1-exp(1)+2;
157  * x1 + x2 - x2*x3+5
158  * ]
159  */
174 void get_non_lin_sys_f3_(vector_t x_arr[], vector_t f3_vec[])
175 {
176  double complex z;
177 
178  f3_vec[0] = exp(-x_arr[0] * x_arr[1]) + log(x_arr[0]) - exp(-2);
179  if (x_arr[2] < 0) {
180  z = csqrtl(x_arr[2]);
181  f3_vec[1] = exp(x_arr[0]) - creal(z) / x_arr[0] - exp(1) + 2;
182  }
183  else {
184  f3_vec[1] = exp(x_arr[0]) - sqrt(x_arr[2]) / x_arr[0] - exp(1)
185  + 2;
186  }
187 
188  f3_vec[2] = x_arr[0] + x_arr[1] - x_arr[1] * x_arr[2] + 5;
189 }
190 
191 /*
192  * [ 1/x1 - x2*exp(-x1*x2), -x1*exp(-x1*x2), 0]
193  * Jf3 = [ exp(x1) + x3^(1/2)/x1^2, 0, -1/(2*x1*x3^(1/2))]
194  * [ 1, 1 - x3, -x2]
195  */
210 void get_non_lin_sys_J3_(vector_t x_arr[], matrix_t J3[][3])
211 {
212  double complex z;
213 
214  J3[0][0] = 1 / x_arr[0] - x_arr[1] * exp(-x_arr[0] * x_arr[1]);
215  J3[0][1] = -x_arr[0] * exp(-x_arr[0] * x_arr[1]);
216  J3[0][2] = 0;
217  if (x_arr[2] < 0) {
218  z = csqrtl(x_arr[2]);
219  J3[1][0] = exp(x_arr[0]) + creal(z) / pow(x_arr[0], 2);
220  }
221  else {
222  J3[1][0] = exp(x_arr[0]) + sqrt(x_arr[2]) / pow(x_arr[0], 2);
223  }
224  J3[1][1] = 0;
225  if (x_arr[2] < 0) {
226  z = csqrtl(x_arr[2]);
227  if ((x_arr[0] == 0) || (creal(z) == 0)) {
228  //J3[1][2] = -1;
229  J3[1][2] = -DBL_MAX;
230  }
231  else {
232  J3[1][2] = -1 / (2 * x_arr[0] * creal(z));
233  }
234 
235  }
236  else {
237  J3[1][2] = -1 / (2 * x_arr[0] * sqrt(x_arr[2]));
238  }
239 
240  J3[2][0] = 1;
241  J3[2][1] = 1 - x_arr[2];
242  J3[2][2] = -x_arr[1];
243 }
244 
245 void fsolve_test(void)
246 {
247  puts("############ Test Solve Non-linear Algebra Module ############");
248  puts(
249  "\n******************************* f1 ******************************");
250  vector_t f1_vec[2];
251  vector_t x1_arr[2] = { 3, 2 };
252  matrix_t J1[2][2];
253  get_non_lin_sys_f1_(x1_arr, f1_vec);
254  get_non_lin_sys_J1_(x1_arr, J1);
255  uint8_t f1_length = 2;
256  vector_t x0_arr1[2] = { .5, .5 };
257  uint8_t n = 2;
258  vector_t est_x1[n];
259  uint8_t iter_num;
260  enum NON_LIN_ALGORITHM algo;
261 
262  algo = Newton_Raphson;
263  iter_num = fsolve(f1_length, n, x0_arr1, algo, est_x1,
265 // iter_num = newton_raphson(f1_length, n, x0_arr1, eps, max_it_num, est_x1,
266 // &get_non_lin_sys_f1_, &get_non_lin_sys_J1_);
267  printf("iter_num = %u\nx1 = ", iter_num);
268  vector_flex_print(n, est_x1, 3, 3);
269  puts("");
270  get_non_lin_sys_f1_(est_x1, f1_vec);
271  printf("f1(x1) = ");
272  vector_flex_print(f1_length, f1_vec, 3, 3);
273 
274  puts("\nDamped Newton");
275  iter_num = 0;
276  vector_clear(n, est_x1);
277  vector_clear(f1_length, f1_vec);
278  algo = Damped_Newton_Raphson;
279  iter_num = fsolve(f1_length, n, x0_arr1, algo, est_x1,
281  printf("iter_num = %u\nx1 = ", iter_num);
282  vector_flex_print(n, est_x1, 3, 3);
283  puts("");
284  get_non_lin_sys_f1_(est_x1, f1_vec);
285  printf("f1(x1) = ");
286  vector_flex_print(f1_length, f1_vec, 3, 3);
287 
288  puts(
289  "/\n******************************** f2 ******************************/");
290  vector_t f2_vec[3];
291  matrix_t J2[3][3];
292  vector_t x2_arr[3] = { 0.1, 0.2, 0.3 };
293  get_non_lin_sys_f2_(x2_arr, f2_vec);
294  get_non_lin_sys_J2_(x2_arr, J2);
295  uint8_t f2_length = 3;
296  n = 3;
297  vector_t est_x2[n];
298  vector_t x0_arr2[3] = { 1, 1, 1 };
299  algo = Newton_Raphson;
300  iter_num = fsolve(f2_length, n, x0_arr2, algo, est_x2,
302 // iter_num = newton_raphson(f2_length, n, x0_arr2, eps, max_it_num, est_x2,
303 // &get_non_lin_sys_f2_, &get_non_lin_sys_J2_);
304  printf("iter_num = %u\nx2 = ", iter_num);
305  vector_flex_print(n, est_x2, 3, 9);
306  puts("");
307  get_non_lin_sys_f2_(est_x2, f2_vec);
308  printf("f2(x2) = ");
309  vector_flex_print(f2_length, f2_vec, 3, 3);
310  puts("");
311 
312  puts("\nDamped Newton");
313  iter_num = 0;
314  vector_clear(n, est_x2);
315  vector_clear(f2_length, f2_vec);
316  algo = Damped_Newton_Raphson;
317  iter_num = fsolve(f2_length, n, x0_arr2, algo, est_x2,
319  printf("iter_num = %u\nx2 = ", iter_num);
320  vector_flex_print(n, est_x2, 3, 9);
321  puts("");
322  get_non_lin_sys_f2_(est_x2, f2_vec);
323  printf("f2(x2) = ");
324  vector_flex_print(f2_length, f2_vec, 3, 3);
325 
326  puts(
327  "/\n******************************** f3 ******************************/");
328 
329  vector_t f3_vec[3];
330  vector_t x3_arr[3] = { 1e1, 2e2, 3e3 };
331  matrix_t J3[3][3];
332  get_non_lin_sys_f3_(x3_arr, f3_vec);
333  get_non_lin_sys_J3_(x3_arr, J3);
334 
335  uint8_t f3_length = 3;
336  n = 3;
337  vector_t est_x3[n];
338  vector_t x0_arr3[3] = { 1, 1, 1 };
339  algo = Newton_Raphson;
340  iter_num = fsolve(f3_length, n, x0_arr3, algo, est_x3,
343 // iter_num = newton_raphson(f3_length, n, x0_arr3, eps, max_it_num, est_x3,
344 // &get_non_lin_sys_f3_, &get_non_lin_sys_J3_);
345  printf("iter_num = %u\nx2 = ", iter_num);
346  vector_flex_print(n, est_x3, 3, 9);
347  puts("");
348  get_non_lin_sys_f3_(est_x3, f3_vec);
349  printf("f3(x3) = ");
350  vector_flex_print(f3_length, f3_vec, 3, 3);
351  puts("");
352 
353  puts("\nDamped Newton");
354  iter_num = 0;
355  vector_clear(n, est_x3);
356  vector_clear(f3_length, f3_vec);
357  iter_num = fsolve(f3_length, n, x0_arr3, algo, est_x3,
360  printf("iter_num = %u\nx3 = ", iter_num);
361  vector_flex_print(n, est_x3, 3, 9);
362  puts("");
363  get_non_lin_sys_f3_(est_x3, f3_vec);
364  printf("f3(x3) = ");
365  vector_flex_print(f3_length, f3_vec, 3, 3);
366  puts("");
367 
368 }
newton_raphson.h
Implement the Newton–Raphson algorithm.
get_non_lin_sys_J1_
void get_non_lin_sys_J1_(vector_t x_arr[], matrix_t J1[][2])
The Jacobian matrix of the non-linear system .
Definition: fsolve_test.c:80
vector_flex_print
void vector_flex_print(uint32_t length, vector_t arr[], uint8_t before_dot, uint8_t after_dot)
Display the values of the vector's elements.
Definition: vector.c:284
fsolve.h
Solve multi-variant nonlinear equation systems.
get_non_lin_sys_J3_
void get_non_lin_sys_J3_(vector_t x_arr[], matrix_t J3[][3])
The Jacobian matrix of the non-linear system .
Definition: fsolve_test.c:210
vector_t
#define vector_t
Define the data type of the vector elements.
Definition: vector.h:33
get_non_lin_sys_f1_
void get_non_lin_sys_f1_(vector_t x_arr[], vector_t f1_vec[])
The non-linear system f1.
Definition: fsolve_test.c:54
NON_LIN_ALGORITHM
NON_LIN_ALGORITHM
Possible algorithms to solve multi-variant nonlinear equation systems.
Definition: fsolve.h:30
Damped_Newton_Raphson
Damped Newton–Raphson algorithm.
Definition: fsolve.h:32
matrix.h
Matrix computations.
damped_newton_raphson.h
Implement the damped Newton–Raphson algorithm.
M_PI
#define M_PI
Definition: matrix.h:54
vector_clear
void vector_clear(uint8_t size, vector_t arr[])
Clear all the elements of the vector.
Definition: vector.c:32
get_non_lin_sys_J2_
void get_non_lin_sys_J2_(vector_t x_arr[], matrix_t J2[][3])
The Jacobian matrix of the non-linear system .
Definition: fsolve_test.c:139
get_non_lin_sys_f2_
void get_non_lin_sys_f2_(vector_t x_arr[], vector_t f2_vec[])
The non-linear system .
Definition: fsolve_test.c:110
fsolve_test
void fsolve_test(void)
Examples of solving non-linear equation systems.
Definition: fsolve_test.c:245
Newton_Raphson
Newton–Raphson algorithm.
Definition: fsolve.h:31
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
get_non_lin_sys_f3_
void get_non_lin_sys_f3_(vector_t x_arr[], vector_t f3_vec[])
The non-linear system .
Definition: fsolve_test.c:174
vector.h
Vector computations.
fsolve
uint8_t fsolve(uint8_t f_length, uint8_t x0_length, vector_t x0_arr[], enum NON_LIN_ALGORITHM algo, 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[][x0_length]))
Solve systems of multi-variant nonlinear equations.
Definition: fsolve.c:29