Численное решение задачи Коши для обыкновенных дифференциальных уравнений методами Рунге-Кутта

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

ПЕНЗЕНСКИЙ  ГОСУДАРСТВЕННЫЙ  УНИВЕРСИТЕТ

ПОЛИТЕХНИЧЕСКИЙ  ИНСТИТУТ

Отчет о выполнении лабораторной работы №3

«Численное решение задачи Коши для обыкновенных дифференциальных уравнений методами Рунге-Кутта»

по дисциплине

«Численные методы»

Выполнил: студент гр. 19ВГ1

Авоськин П.И

Проверила: Кудряшова Н.Ю.

Пенза, 2021

Цель работы. Решение обыкновенного дифференциального уравнения первого порядка методами Рунге-Кутта.

Задание: написать программу для решения обыкновенногодифференциального  уравнения первого порядка методами Рунге-Кутта.

Варианты заданий

[pic 1]

Теоретическая часть

Пусть требуется найти на отрезке [x0,x0+X] решение дифференциального уравнения

[pic 2]

Разобьем [x0,x0+X] на отрезки [xj-1,xj]. Последовательно будем получать приближения yj к значениям решения yj(x), j = 1,..., N. Пусть значение yj уже найдено, тогда значение yj+1 будем определять по следующей расчетной формуле

[pic 3]

Этот метод называется методом Эйлера. Можно построить другой класс расчетных формул, к которому принадлежит метод Эйлера. Укажем сначала простейшие методы этого класса, получаемые из наглядных соображений. Пусть известно значение решения y(x) и требуется вычислить значение y(x+h). Рассмотрим равенство

[pic 4]

При замене интеграла в правой части на величинупогрешность [pic 5]

имеет порядок , т.е.  поскольку  то отсюда имеем  Отбрасывая член порядка  и обозначая xj = x, x + h = xj+1, получим расчетную формулу Эйлера (2). Для получения более точной расчетной формулы нужно точнее аппроксимировать интеграл в правой части (3). Воспользовавшись формулой трапеций, получим[pic 6][pic 7][pic 8][pic 9][pic 10]

,[pic 11]

иначе

[pic 12]

Соответствующая расчетная формула:

[pic 13]

Обычно это уравнение неразрешимо явно относительно yj+1. Поэтому на

практике чаще всего пользуются расчетными формулами

[pic 14]

Если интеграл в правой части (3) заменить по формуле прямоугольников, то можно построить другую пару расчетных формул с погрешностью на шаге того же порядка:

[pic 15]

Полученные методы относятся к семейству методов Рунге-Кутта.

Листинг программы:

#include<iostream>

#include<iomanip>

#include<cmath>

using namespace std;

double a = 0.1, b = 2, y_0 = 1;

int n = 20;

long double foo(double x, double y) {

    return y/x;

    //return atan(x)/x + y;

}

void first() {

    double h = (b - a) / n;

    double h1 = h / 2;

    long double* x = new long double[n + 1];

    for (int i = 0; i <= n; i++)

    {

        x[i] = a + h * i;

    }

    long double* x1 = new long double[2 * n + 1];

    for (int i = 0; i <= 2 * n; i++)

    {

        x1[i] = a + h1 * i;

    }

    long double* y = new long double[n + 1];

    y[0] = y_0;

    for (int i = 0; i <= n - 1; i++)

    {

        y[i + 1] = y[i] + h * foo(x[i], y[i]);