Lab12/Lab_12/Source.c - master - aytsan_ns/Num_meth

Lab12/Lab_12/Source.c
190 строк4 KB
aytsan-ns
Изменения в код 12 работы. Убрала повторяющееся вычисление k1, изменила отчет
30 сен 2024, 18:38
30 сен 2024, 18:383794d15
#define _CRT_SECURE_NO_WARNINGS
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#define PI 3.14159265358979323846

double function(double x, double y) {
	return (y / x + x * cos(x));
}

double F(double x) {
	return x * sin(x);
}

double RungeKutta(double x0, double y0, double h) {
	double y = y0;
	double x = x0;

	double k1 = function(x, y);
	double k2 = function(x + h / 2, y + (h / 2) * k1);
	double k3 = function(x + h, y - h * k1 + 2 * h * k2);

	y = y + (h / 6) * (k1 + 4 * k2 + k3);
	return y;
}

double* RungeKuttaEps(double x0, double y0, double xn, double eps, int* n, double** error_rate, double** result_h) {

	double h = xn - x0;
	double* result = NULL;
	double x1, x2, y1, y2_0, y2, k1, k2, k3;

	while (x0 < xn) {
		if (x0 + h > xn) h = xn - x0;

		x1 = x0;
		y1 = y0;

		k1 = function(x0, y0);

		k2 = function(x1 + h / 2, y1 + (h / 2) * k1);
		k3 = function(x1 + h, y1 - h * k1 + 2 * h * k2);
		y1 = y1 + (h / 6) * (k1 + 4 * k2 + k3);
		x1 = x1 + h;

		do {
			x2 = x0;
			y2_0 = y0;

			k2 = function(x2 + h / 4, y2_0 + (h / 4) * k1);
			k3 = function(x2 + h / 2, y2_0 - (h / 2) * k1 + h * k2);
			y2_0 = y2_0 + (h / 12) * (k1 + 4 * k2 + k3);
			x2 = x2 + h / 2;

			k1 = function(x2, y2_0);
			k2 = function(x2 + h / 4, y2_0 + (h / 4) * k1);
			k3 = function(x2 + h / 2, y2_0 - (h / 2) * k1 + h * k2);
			y2 = y2_0 + (h / 12) * (k1 + 4 * k2 + k3);
			x2 = x2 + h / 2;

			double actual_error = fabs(y2 - y1) / 7;

			if (actual_error < eps) {
				*n += 1;

				double* new_error = realloc((*error_rate), sizeof(double) * (*n));
				if (new_error == NULL) {
					printf("Error");
					return NULL;
				}
				*error_rate = new_error;
				(*error_rate)[(*n) - 1] = actual_error;

				double* new_result = realloc(result, sizeof(double) * (*n));
				if (new_result == NULL) {
					printf("Error");
					return NULL;
				}
				result = new_result;
				result[(*n) - 1] = y2;

				double* new_result_h = realloc((*result_h), sizeof(double) * (*n));
				if (new_result_h == NULL) {
					printf("Error");
					return NULL;
				}
				*result_h = new_result_h;
				(*result_h)[(*n) - 1] = h / 2;

				x0 = x2;
				y0 = y2;

				if (actual_error < eps / 2) h *= 2;

				break;
			}
			else {
				y1 = y2_0;
				h /= 2;
			}
		} while (1);
	}
	return result;
}

double MaxError(double* error, int n) {
	double max_error = 0;
	for (int i = 0; i < n; i++)
		if (error[i] > max_error) max_error = error[i];
	return max_error;
}

void Task1() {
	double a = PI / 2, b = 2 * PI, xi = PI / 2, yi = PI / 2;
	int n1 = 3, n2 = 9;
	double h1 = (b - a) / n1, h2 = (b - a) / n2;
	FILE* result_h1 = fopen("result_h1.txt", "w");
	FILE* result_h2 = fopen("result_h2.txt", "w");
	FILE* error_result_h1 = fopen("error_result_h1.txt", "w");
	FILE* error_result_h2 = fopen("error_result_h2.txt", "w");

	fprintf(result_h1, "%.10lf %.10lf\n", xi, yi);
	fprintf(error_result_h1, "%.10lf %.10lf\n", xi, fabs(yi - F(xi)));
	for (int i = 1; i <= n1; i++) {
		yi = RungeKutta(xi, yi, h1);
		xi = a + h1 * i;
		fprintf(result_h1, "%.10lf %.10lf\n", xi, yi);
		fprintf(error_result_h1, "%.10lf %.10lf\n", xi, fabs(yi - F(xi)));
	}

	xi = PI / 2;
	yi = PI / 2;
	fprintf(result_h2, "%.10lf %.10lf\n", xi, yi);
	fprintf(error_result_h2, "%.10lf %.10lf\n", xi, fabs(yi - F(xi)));
	for (int i = 1; i <= n2; i++) {
		yi = RungeKutta(xi, yi, h2);
		xi = a + h2 * i;
		fprintf(result_h2, "%.10lf %.10lf\n", xi, yi);
		fprintf(error_result_h2, "%.10lf %.10lf\n", xi, fabs(yi - F(xi)));
	}

	_fcloseall();
}

void Task2() {
	double a = PI / 2, b = 2 * PI;
	double x = a, y = PI / 2;
	double eps = 1e-7;
	FILE* result_h = fopen("result_h.txt", "w");

	int n = 0;
	double* error = NULL;
	double* h = NULL;
	double* result = RungeKuttaEps(a, y, b, eps, &n, &error, &h);

	double H = 0;

	for (int i = 0; i < n; i++) {
		H += (h[i] * 2);
		fprintf(result_h, "%.10lf %.10lf\n", x, h[i]);
		x = x + 2 * h[i];
		fprintf(result_h, "%.10lf %.10lf\n", x, h[i]);
	}

	_fcloseall();
}

void Task3() {
	double a = PI / 2, b = 2 * PI;
	double x = a, y = PI / 2;
	FILE* error_eps = fopen("error_eps.txt", "w");

	for (double eps = 0.1; eps >= 1e-7; eps /= 10) {
		int n = 0;
		double* error = NULL;
		double* h = NULL;
		double* result = RungeKuttaEps(a, y, b, eps, &n, &error, &h);
		double max_error = MaxError(error, n);
		fprintf(error_eps, "%.10lf %.10lf\n", eps, max_error);
	}

	_fcloseall();
}

void main() {
	Task1();
	Task2();
	Task3();
}