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

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Санкт-Петербургский государственный политехнический университет ИПММ

Кафедра Телематики

Лабораторная работа

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

Работу выполнил:

Подорова А.Я.

группа 23607/1

Работу проверил:

Кадырова Н.О.



Цель работы

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

Данное уравнение:

Результатом решения данного уравнения получено следующее точное значение: .

Решить краевую задачу для дифференциального уравнения второго порядка путем сведения краевой задачи к задаче Коши

Данное уравнение:

Результатом решения данного уравнения получено следующее точное значение:

Метод Эйлера-Коши (Рунге-Кутты второго порядка) для ДУ первого порядка.

Численные методы решения задачи Коши

,

на равномерной сетке отрезка с шагом являются методами Рунге-Кутта, если, начиная с данных решение ведётся по следующим рекуррентным формулам:



(8.6)

Метод называют методом Рунге-Кутта порядка если он имеет -й порядок точности по шагу на сетке. Порядок точности достигается с помощью формул (8.6) при определённых значениях коэффициентов и коэффициент всегда полагают равным нулю. Эти коэффициенты вычисляют по следующей схеме:

1) точное решение и его приближение представляют в виде разложения по формуле Тейлора с центром вплоть до слагаемого порядка

2) из равенств подобных членов при одинаковых степенях в двух разложениях получают уравнения, решая которые находят коэффициенты и .

Метод Рунге-Кутта второго порядка называют методом Эйлера-Коши, если Алгоритм метода Эйлера-Коши следует из формул (8.6):

(8.7)

Для практической оценки погрешности решения можно применять правило Рунге, полагая в формуле (8.5) 

Решение краевой задачи для ДУ второго порядка.

Дано линейное дифференциальное уравнение второго порядка:

(1.1)

с краевыми условиями:

(1.2)

где функции P(x), Q(x), F(x) непрерывны, и - заданные постоянные, причем: и .

Требуется найти решение y(x), удовлетворяющее краевым условиям (1.2).

Описание метода

Решение дифференциального уравнения (1.1) с краевыми условиями (1.2) будем искать в виде линейной комбинации:

(2.1)

где c=const

Подставим y(x) в виде (2.1) в исходное дифференциальное уравнение (1.1) и получим выражение:

(2.2)

Необходимо, чтобы коэффициенты при постоянной с обращались в ноль, получим систему уравнений:



(2.3)

Из (2.3) видно, что функция u=u(x) - ненулевое решение соответствующего однородного уравнения, а v=v(x) - некоторое решение данного неоднородного уравнения (1.1).

Чтобы свести краевую задачу к задачам Коши для функций u=u(x) и v=v(x), подставим в первое краевое условие (1.2) выражение для функции y(x) и получим:

(2.4)

Для того чтобы равенство (2.4) было справедливо при любом с, необходимо и достаточно, чтобы коэффициенты при постоянной с обращались в ноль, т. е. должны быть выполнены равенства:

(2.5)

Получили систему (2.5) из двух уравнений с четырьмя неизвестными. Решений системы будет бесконечное множество. Найдем хотя бы одно.

Для обеспечения первого равенства системы, например, можно подобрать:

(2.6),

где постоянная k - отлична от нуля, так как тривиальное решение u(a)=0 можно отбросить.

Для выполнения второго равенства системы (2.5) можно положить

(2.7)

Подберем теперь постоянную c так, чтобы функция y(x) удовлетворяла краевому условию (1.2) на конце x=b. Это дает:

откуда:

,

при этом предполагается, что знаменатель

Для решения полученных уравнений из системы (2.3) будем использовать метод Рунге-Кутта, который имеет достаточно высокую точность на всем интервале порядка

Рассмотрим второе дифференциальное уравнение из системы (2.3) с начальным условием (2.7). Дифференциальные уравнения системы (2.3) являются однотипными, поэтому решение первого уравнения системы (2.3) с начальным условием (2.6) осуществляется аналогичным способом, при условии .

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

и подставим ее во второе дифференциальное уравнение из системы (2.3), получим:



,

Расчет производим следующим образом: выберем шаг h, приращение x в зависимости от шага будет:

где n=0,1…k-1

Соответствующие значения и искомых функций v и z определяются формулами:

Метод Рунге-Кутты для ДУ первого порядка

Результат работы программы для заданных 10 разбиений.



x	RK2	ansDY	eps
2	4	4	0
2,4	5,7611	5,76	0,0011
2,8	7,8387	7,84	0,0013
3,2	10,2342	10,24	0,0058
3,6	12,9484	12,96	0,0116
4	15,9816	16	0,0184
4,4	19,3338	19,36	0,0262
4,8	23,0053	23,04	0,0347
5,2	26,996	27,04	0,044
5,6	31,3061	31,36	0,0539
6	35,9354	36	0,0646

Получены следующие результаты:

В первом столбце приведено разбиение отрезка [2;6] с шагом h=(b-a)/n , что для n=10 есть 0.4;

Во стором столбце значение функции, вычисленное по методу Эйлера-Коши.

В третьем столбце приведено точное значение, соответствующее формуле y=x², что является решением дифференциального уравнения.

В четвертом столбце приведена погрешность проведенных вычисленний.(eps=|yi-y(x)|).

Опытным путем было установлено, что при увеличении числа разбиений n погрешность метода Эйлера-Коши минимизируется.

Можно подтвердить, привидя пример для n=20, и проанализировать:

x	RK2	ansDY	eps
2	4	4	0
2,2	4.840495867768595	4.840	4.958677685946711e-004
2,4	5.760626930461642	5.760	6.269304616406402e-004
2,6	6.760508252367625	6.760	5.082523676227524e-004
2,8	7.840207546326258	7.840	2.075463262540822e-004
3	8.999767031613843	9.00	2.329683861628240e-004
3,2	10.239214312181240	10.240	7.856878187677552e-004
3,4	11.558568164571458	11.560	0.001431835428550
3,6	12.957841794432600	12.960	0.002158205567410
3,8	14.437044757653958	14.44	0.002955242346054
4	15.996184140722654	16.00	0.003815859277360
4,2	17.635265311934894	17.64	0.004734688065120
4,4	19.354292414432642	19.36	0.005707585567375
4,6	21.153268698687814	21.16	0.006731301312207
4,8	23.032196752172062	23.040	0.007803247827962
5	24.991078661444632	25.00	0.008921338555396
5,2	27.029916128766249	27.040	0.010083871233782
5,4	29.148710557463726	29.16	0.011289442536306
5,6	31.347463115407390	31.36	0.012536884592645
5,8	33.626174782891027	33.64	0.013825217109009
6	35.984846389219804	36.00	0.015153610780239

Сравнивая привиденные таблицы, делаем вывод, что при увеличении n погрешность метода уменьшается.

Решение краевой задачи для ДУ второго порядка

Результат работы программы для 10 и 20 разбиений:



В первом столбце нумерация узла.

Во втором столбце приведено разбиение отрезка [0; pi/2] с шагом h = 0.16, в случае 10 разбиений и с шагом h = 0.08 в случае 20 разбиений.

В третьем столбце значение функции y(x_n), полученное запрограммированным методом.

В четвертом столбце значение функции sin(x_n), которое является решением данного ДУ.

В последнем столбце значения найденных абсолютных погрешностей.

Аналогично методу Эйлера-Коши, при увеличении числа разбиений n, погрешность уменьшается.

Вывод:

При проведении данной работы подтвердилось, что точность решения и в случае решения краевой задачи для ДУ второго порядка, и при применении метода Рунге-Кутты для ДУ первого порядка зависит от количества разбиений. Чем больше количество разбиений, тем точность выше.

Код программ

Решения краевой задачи для ДУ второго порядка

#include "iostream"

#include "math.h"

#include <locale.h>

#include <conio.h>

using namespace std;

const int Nmax = 100;

# define M_PI_2 1.57079632679489661923 /* pi/2 */

const double a = 0, b = M_PI_2, A = 0, B = 1, epsel = 0.00001;

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

double Q(double x){ return (3.);}

double P(double x){ return -tan(x);}

void RK4(double h, int n, double U[], double V[], double z[], double t[], double a)

{

double k1z, k2z, k3z, k4z,

k1v, k2v, k3v, k4v,

k1t, k2t, k3t, k4t,

k1u, k2u, k3u, k4u;

double x = a;

z[0] = 0;

V[0] = A;

t[0] = 1;

U[0] = 0;

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

{

x = a + i*h;

k1v = z[i];

k1z = F(x) - P(x)*z[i] - Q(x)*V[i];

k2v = z[i] + h/2 * k1z;

k2z = F(x + h/2) - P(x + h/2)*(z[i] + h/2 * k1z) - Q(x + h/2)*(V[i] + h/ 2 * k1v);

k3v = z[i] + h / 2 * k2z;

k3z = F(x + h / 2) - P(x + h / 2)*(z[i] + h / 2 * k2z) - Q(x + h / 2)*(V[i] + h / 2 * k2v);

k4v = z[i] + h*k3z;

k4z = F(x + h) - P(x + h)*(z[i] + h*k3z) - Q(x + h)*(V[i] + h*k3v);

k1u = t[i];

k1t = -P(x)*t[i] - Q(x)*U[i];

k2u = t[i] + h / 2 * k1t;

k2t = -P(x + h / 2)*(t[i] + h / 2 * k1t) - Q(x + h / 2)*(U[i] + h / 2 * k1u);

k3u = t[i] + h / 2 * k2t;

k3t = -P(x + h / 2)*(t[i] + h / 2 * k2t) - Q(x + h / 2)*(U[i] + h / 2 * k2u);

k4u = t[i] + h*k3t;

k4t = -P(x + h)*(t[i] + h*k3t) - Q(x + h)*(U[i] + h*k3u);

V[i + 1] = V[i] + h / 6 * (k1v + 2 * k2v + 2 * k3v + k4v);

z[i + 1] = z[i] + h / 6 * (k1z + 2 * k2z + 2 * k3z + k4z);

U[i + 1] = U[i] + h / 6 * (k1u + 2 * k2u + 2 * k3u + k4u);

t[i + 1] = t[i] + h / 6 * (k1t + 2 * k2t + 2 * k3t + k4t);

}

}

int main()

{

while (true)

{

setlocale(LC_ALL, "Rus");

int i = 0, j = 0, n = 0, N;

cout << "\nВведите количество разбиений N: ";

cin>>N;

double U[Nmax], V[Nmax], z[Nmax], y[Nmax], yt[Nmax], yv[Nmax];

double h = 0, x = 0, delta = 0, t[Nmax], eps[Nmax], C = 0;

n = N;

do

{

if (n >= Nmax)

{

cout <<"Выход за пределы массива";

return 0;

}

h = (b-a)/n;

RK4(h, n, U, V, z, t, a);

C = (B-V[n])/U[n];

//Нахождение у

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

y[i] = C*U[i] + V[i];

//Проверка точности

delta = 0;

for (i = 0; i <= n / 2; i++)

{

eps[i] = fabs(y[2 * i] - yt[i]);

if (eps[i]>delta)

delta = eps[i];

}

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

{

yv[i] = yt[i];

yt[i] = y[i];

}

n = 2 * n;

}

while (delta>epsel);

j = 0;

x = a;

n=n/4;

int k = n / N;

h = (b - a)/ n;

printf("\tx\ty\t\tsin(x)\t\tAbsEps\n");

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

{

printf("%d\t%1.2f\t%f\t%f\t%1.2e\n", j++, x, yv[i], F(x), eps[i]);

x = x + k*h;

}

cout <<"Для выхода нажмите клавишу 'Esc', для продолжения любую другую клавишу\n";

char yes=_getch();

if(yes==27) return 0;

Метод Рунге-Кутты для ДУ первого порядка

function f1=df(x,y)

f1=5 + (2*x - 5)/x^2*y;

end

function ft=dft(x)

ft=x^2;

end

a=2; b=6; yo=4; n=20;

fu=@df;

fut=@dft;

h=(b-a)/n;

x=a; y=yo;

emax=0; %макс. отклонение

for i=1:n

z=y+h*fu(x,y);

y=y+0.5*h*(fu(x,y)+fu(x+h,z))

x=x+h

eps=abs(y-fut(x))

fut(x)

if eps>emax

emax=eps;

end

дифференциальный уравнение рунге коши

Размещено на Allbest.ru

...