kursach_for_dosha

// diff_solve.cpp : Этот файл содержит функцию "main". Здесь начинается и заканчивается выполнение программы.
//

#include "pch.h"
#include <iostream>
#include <cstdlib>
#pragma warning(disable : 4996)

//Заполнение массива функций
float* FuncSet(float t, float *x, int Nequations)
{
    float *diff_equation;
    diff_equation = (float *)malloc(Nequations * sizeof(float));
    diff_equation[0] = x[0];
    diff_equation[1] = x[1];
    diff_equation[2] = x[2];
    return diff_equation;
}

//Функция вывода на экран и печати в файл
void Printing(float t, float *y, int N_eqs, FILE* output)
{
    fprintf(output, "%2.2f\t", t);
    printf("t=%2.2f\t", t);
    for (int i = 0; i < N_eqs; i++)
    {
        fprintf(output, "%2.2f\t", y[i]);
        printf("y%i=%2.2f\t", i, y[i]);
    }
    fprintf(output, "\n");
    printf("\n");
}

//Вычисление методом Эйлера на одном шаге для каждой функции
float *StepE(float *y, float t, float step, float* functions(float, float*, int), int N_eqs)
{
    for (int k = 0; k < N_eqs; k++)
    {
        y[k] = y[k] + step * functions(t, y, N_eqs)[k];
    }
    return y;
}

//Вычисление методом Рунге-Кутта на одном шаге для каждой функции
float *StepRK(float *y, float t, float step, float* functions(float, float*, int), int N_eqs)
{
    //Объявление К0, К1, К2 и К3
    float *K0 = (float*)malloc(N_eqs * sizeof(float));
    float *K1 = (float*)malloc(N_eqs * sizeof(float));
    float *K2 = (float*)malloc(N_eqs * sizeof(float));
    float *K3 = (float*)malloc(N_eqs * sizeof(float));

    //создание массива значений n-1 шага
    float* step_meaning = (float*)malloc(N_eqs * sizeof(float));

    for (int i = 0; i < N_eqs; i++)
    {
        step_meaning[i] = y[i];
    }

    //Объявление вспомогательного массива
    float *current_RK = (float*)malloc((N_eqs) * sizeof(float));

    int j;

    //Циклы расчёта значений вспомогательных аргументов
    //Вычисление всех К0 и у+К0/2

        for (j = 0; j < N_eqs; j++)
        {
            K0[j] = step * functions(t, step_meaning, N_eqs)[j];
            current_RK[j] = step_meaning[j] + K0[j] / 2;
        }

        //Вычисление всех К1 и у+К1/2
        for (j = 0; j < N_eqs; j++)
        {
            K1[j] = step * functions((t + step / 2), (current_RK), N_eqs)[j];
            current_RK[j] = step_meaning[j] + K1[j] / 2;
        }

        //Вычисление К2 и у+К2
        for (j = 0; j < N_eqs; j++)
        {
            K2[j] = step * functions((t + step / 2), (current_RK), N_eqs)[j];
            current_RK[j] = step_meaning[j] + K2[j];
        }

        //Вычисление К3
        for (j = 0; j < N_eqs; j++)
        {
            K3[j] = step * functions((t + step), (current_RK), N_eqs)[j];
        }

    //Цикл расчёта знaчений функций
    for (j = 0; j < N_eqs; j++)
    {
        y[j] = y[j] + (K0[j] + 2 * K1[j] + 2 * K2[j] + K3[j]) / 6;
    }

    return y;
}

//Возврат массива к начальным условиям
void recall(float*y, float* y0, int N_eqs)
{
    for (int i = 0; i < N_eqs; i++)
    {
        y[i] = y0[i];
    }
}


int main()
{
    int N_eqs = 3;

    //Объявление массива функций
    float (**diff_equation)(float,float*) = NULL;
    diff_equation=(float(**)(float, float*))malloc((N_eqs) * sizeof(float*(*)(float, float*)));

    //объявление массива решений на шаге
    //[0] - значение у1
    //[1] - значение у2
    //[2] - значение у3
    float *y = (float*)malloc((N_eqs) * sizeof(float));

    //Начальное заполнение массива решений
    y[0] = 1;
    y[1] = 1;
    y[2] = 1;

    //Создание массива начальных значений
    float* y0 = (float*)malloc((N_eqs) * sizeof(float));
    for (int i=0; i < N_eqs; i++)
    {
        y0[i] = y[i];
    }

    //t - переменная
    //step - шаг системы
    //а - меншая граница диапазона
    //b - большая граница диапазона
    float t , a = 0, b = 5, step = 0.1F;

    //Объявление файлов вывода
    FILE* output_Eiler;
    FILE* output_RK;

    t = a;

    //Открытие/создание файлов
    output_Eiler = fopen("output_Eiler.txt", "w");
    output_RK = fopen("output_RK.txt", "w");

    //Вызов функции заполнения массива значений

    printf("Calculated by Eiler method.\n");

    Printing(t, y, N_eqs, output_Eiler);
    while (t < (b - step))
    {
        t = t + step;
        y = StepE(y, t, step, FuncSet, N_eqs);
        Printing(t, y, N_eqs, output_Eiler);
    }

    t = a;

    recall(y, y0, N_eqs);

    printf("Calculated by RK method.\n");

    Printing(t, y, N_eqs, output_RK);
    while (t < (b - step))
    {
        t = t + step;
        y = StepRK(y, t, step, FuncSet, N_eqs);
        Printing(t, y, N_eqs, output_RK);
    }

    free(y);
    free(y0);

    fclose(output_Eiler);
    fclose(output_RK);
}