NTNU-Numerical-Method/wk6-nonlinear.c at main · KennethOng02/NTNU-Numerical-Method · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "matrix.h"

typedef struct {
  double h;
  double k;
  double r;
}jacobian_matrix;

double **getdmatrix1(jacobian_matrix x) {
  double **result = malloc(3 * sizeof(double));
  for(int i = 0; i < 3; i++)
    result[i] = malloc(3 * sizeof(double));
  result[0][0] = 16 + 2 * x.h;
  result[0][1] = 8 + 2 * x.k;
  result[0][2] = -2 * x.r;

  result[1][0] = -12 + 2 * x.h;
  result[1][1] = -18 + 2 * x.k;
  result[1][2] = -2 * x.r;

  result[2][0] = -8 + 2 * x.h;
  result[2][1] = 18 + 2 * x.k;
  result[2][2] = -2 * x.r;

  return result;
}

double **getdmatrix2(jacobian_matrix x) {
  double **result = malloc(3 * sizeof(double));
  for(int i = 0; i < 3; i++)
    result[i] = malloc(3 * sizeof(double));
  result[0][0] = 2 + 2 * x.h;
  result[0][1] = -12 + 2 * x.k;
  result[0][2] = -2 * x.r;

  result[1][0] = 4 + 2 * x.h;
  result[1][1] = 12 + 2 * x.k;
  result[1][2] = -2 * x.r;

  result[2][0] = -10 + 2 * x.h;
  result[2][1] = 2 * x.k;
  result[2][2] = -2 * x.r;

  return result;
}

double *getfmatrix1(jacobian_matrix x) {
  double *result = malloc(3 * sizeof(double));

  result[0] = -1 * (pow((-8 - x.h), 2) + pow((-4 - x.k), 2) - pow(x.r, 2));
  result[1] = -1 * (pow((6 - x.h), 2) + pow((9 - x.k), 2) - pow(x.r, 2));
  result[2] = -1 * (pow((4 - x.h), 2) + pow((-9 - x.k), 2) - pow(x.r, 2));

  return result;
}

double *getfmatrix2(jacobian_matrix x) {
  double *result = malloc(3 * sizeof(double));

  result[0] = -1 * (pow((-1 - x.h), 2) + pow((6 - x.k), 2) - pow(x.r, 2));
  result[1] = -1 * (pow((-2 - x.h), 2) + pow((-6 - x.k), 2) - pow(x.r, 2));
  result[2] = -1 * (pow((5 - x.h), 2) + pow((0 - x.k), 2) - pow(x.r, 2));

  return result;
}

double mdiff(jacobian_matrix initial, jacobian_matrix new) {
  double total = 0;
  total += fabs(initial.h - new.h);
  total += fabs(initial.k - new.k);
  total += fabs(initial.r - new.r);
  return total;
}

int main() {
  int row = 3; int col = 3;

  for(int k = 0; k < 2; k++) {
    jacobian_matrix initial_vector = {2, 2, 2};
    jacobian_matrix new_vector = {1, 1, 1};

    printf("Test %d\n", k + 1);

    while(mdiff(initial_vector, new_vector) >= 10e-6) {
      initial_vector.h = new_vector.h;
      initial_vector.k = new_vector.k;
      initial_vector.r = new_vector.r;

      double **matrix;
      double *oldb;
      if(k == 0) {
        matrix = getdmatrix1(initial_vector);
        oldb   = getfmatrix1(initial_vector);
      }else {
        matrix = getdmatrix2(initial_vector);
        oldb   = getfmatrix2(initial_vector);
      }

      double *x = gaussian_solve(matrix, oldb, row, col);

      new_vector.h = initial_vector.h + x[0];
      new_vector.k = initial_vector.k + x[1];
      new_vector.r = initial_vector.r + x[2];
    }

    printf("h = %.4lf\n", initial_vector.h);
    printf("k = %.4lf\n", initial_vector.k);
    printf("r = %.4lf\n", initial_vector.r);
  }
  return 0;
}