num_methods_sem2/lab13/method.c

#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;
}