123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203 |
- #include <math.h>
- #include <stdlib.h>
- #include <stdio.h>
- #include "header.h"
- double* create_grid_even(int n, double a, double b) {
- double* h = malloc(sizeof(double) * n);
- if (h == NULL) {
- exit(16);
- }
- for (int i = 0; i < n;i++) {
- h[i] = a + ((b - a) / (n - 1)) * i;
- }
- return h;
- }
- double solve_eqution(double* right_term) {
- //слева полагаем, что стоит только неизвестная с коэффициентом 1
- // right_term - массив из двух элементов, на первой позиции - коэф. при неизветсной, на второй - нулевой коэф.
- return right_term[1] / (1 - right_term[0]);
- }
- double* adams_moulten(double* start_points, double a, double b, double(*f)(double, double), int n) {
- double* yy = malloc(sizeof(double) * n);
- double h = (b - a) / (n - 1);
- if (yy == NULL) {
- exit(16);
- }
- for (int o = 0;o < 2;o++) {
- yy[o] = start_points[o];
- }
- for (int i = 2;i < n;i++) {
- double* term = f_eq_form(a + i * h);
- term[0] *= 5 * h / 12;
- term[1] *= 5 * h / 12;
- term[1] += yy[i - 1] + (h / 12) * (8 * f(a + (i - 1) * h, yy[i - 1]) - f(a + (i - 2) * h, yy[i - 2]));
- yy[i] = solve_eqution(term);
- }
- return yy;
- }
- double** runge_adams(double eps, double y0, double a, double b, double(*f)(double, double), int* n_, int* len_h) {
- double* x = malloc(sizeof(double)*2);
- double* y = malloc(sizeof(double)*2);
- if (x == NULL) {
- exit(16);
- }
- if (y == NULL) {
- exit(16);
- }
- long double* h_all = malloc(sizeof(long double));// регистрирует все значения шагов
- long double* h_made = malloc(sizeof(long double)*2);// регистрирует значения шага, соответствующие сделаным шагам
- if (h_all == NULL) {
- exit(16);
- }
- if (h_made == NULL) {
- exit(16);
- }
- double* x2 = malloc(sizeof(double) * 4);
- double* y2 = malloc(sizeof(double) * 4);
- if (x2 == NULL) {
- exit(16);
- }
- if (y2 == NULL) {
- exit(16);
- }
- int k = 1, i = 2;
- double dif, h = (b - a) / 2;
- double y_1 = y0 + h * f(a, y0);
- double y_12 = y0 + (h / 2) * f(a, y0);
- double y_2 = y_12 + (h / 2) * f(a + h / 2, y_12);
- while (fabs(y_2 - y_1) > eps / 100) {
- h /= 2;
- y_1 = y0 + h * f(a + h, y0);
- y_12 = y0 + (h / 2) * f(a + h, y0);
- y_2 = y_12 + (h / 2) * f(a + h / 2, y_12);
- }
- y[0] = y0;
- y[1] = y_1;
- y2[0] = y_12;
- y2[1] = y_2;
- x[0] = a;
- x[1] = a + h;
- x2[0] = a+ h/2;
- x2[1] = a + h;
- h_made[0] = 1;
- h_made[1] = h;
- h = (b - a) / 2;
- h_all[0] = h;
- printf("eps = %e\n", eps);
- double* term;
- while (x[i - 1] + h <= b) {
- h_all = realloc(h_all, sizeof(long double) * (k + 1));
- if (h_all == NULL) {
- exit(16);
- }
- h_all[k] = h;
- term = f_eq_form(x[i - 1] + h);
- term[0] *= 5 * h / 12;
- term[1] *= 5 * h / 12;
- term[1] += y[i - 1] + (h / 12) * (8 * f(x[i - 1], y[i - 1]) - f(x[i - 2], y[i - 2]));
- y_1 = solve_eqution(term);
- term = f_eq_form(x[i - 1] + h / 2);
- term[0] *= 5 * h / 24;
- term[1] *= 5 * h / 24;
- term[1] += y[i - 1] + (h / 24) * (8 * f(x2[1], y2[1]) - f(x2[0], y2[0]));
- y_12 = solve_eqution(term);
- term = f_eq_form(x[i - 1] + h);
- term[0] *= 5 * h / 24;
- term[1] *= 5 * h / 24;
- term[1] += y_12 + (h / 24) * (8 * f(x[i - 1] + h / 2, y_12) - f(x2[ 1], y2[1]));
- y_2 = solve_eqution(term);
- dif = fabs(y_2 - y_1) / 7;
- x2[3] = x[i - 1] + h;
- x2[2] = x[i - 1] + h / 2;
- y2[3] = y_2;
- y2[2] = y_12;
- if (dif > eps/8) {
- h = h / 2;
- }
- else {
- x = realloc(x, sizeof(double) * (i + 1));
- if (x == NULL) {
- exit(16);
- }
- x[i] = x[i - 1] + h;
- y = realloc(y, sizeof(double) * (i + 1));
- if (y == NULL) {
- exit(16);
- }
- y[i] = y_1;
-
- x2[0] = x2[2];
- x2[1] = x2[3];
- y2[0] = y2[2];
- y2[1] = y2[3];
- h_made = realloc(h_made, sizeof(long double) * (i + 1));
- if (h_made == NULL) {
- exit(16);
- }
- h_made[i] = h;
- i += 1;
- }
- k++;
- }
- h_all[k - 1] = h;
- *len_h = k;
- *n_ = i;
- double** ans = malloc(sizeof(double*) * 4);
- ans[0] = x;
- ans[1] = y;
- ans[2] = h_all;
- ans[3] = h_made;
- free(x2);
- free(y2);
- return ans;
- }
- double uglify_initial(double initial, int max_procent, int min_procent) {
- double procent = min_procent + (double)(max_procent - min_procent) * ((double)rand() / (RAND_MAX));
- return initial - procent * initial / 100;
- }
|