Лабораторная работа по численным методам. Лабораторная работа №1. Дисциплина Теория разностных схем Отчет по лабораторной работе 1

Название	Дисциплина Теория разностных схем Отчет по лабораторной работе 1
Анкор	Лабораторная работа по численным методам
Дата	21.03.2020
Размер	274.3 Kb.
Формат файла
Имя файла	Лабораторная работа №1.docx
Тип	Решение #112751
страница	6 из 6

1 2 3 4 5 6

II. Краевая задача для обыкновенного дифференциального уравнения второго прядка

Решается следующая краевая задача для неоднородного ОДУ второго

порядка:

Точное решение:

Задача 5. Конечно-разностный метод

Погрешность в зависимости от числа узлов

Число узлов	Шаг	Погрешность
32	0.03125	0.007149
64	0.015625	0.001751
128	0.0078125	0.000467
512	0,001953125	0.000078

График погрешности

Задача 6. Метод стрельбы

Нелинейное уравнение решается метод секущих. Начальные приближения – границы интервала.

Погрешность в зависимости от числа узлов

Число узлов	Шаг	Погрешность
32	0.03125	0.974050
64	0.015625	0.493350
128	0.0078125	0.248250
512	0,001953125	0.062350

Графики решений

График погрешности

Заключение

В результате проделанной лабораторной работы был изучен теоретический материал необходимый для решения начальных и краевых задач для обыкновенных дифференциальных уравнений.

Для каждой поставленной задачи написана вычислительная программа на языке программирования С++, выполняющая необходимые построения и расчеты.

Приложение А

Метод Эйлера

std::vector> DSolve::Euler(const double N1, const double dt, std::function F, std::function G, double x0, double y0)

{

const size_t N = N1;

std::vector> X;

X.resize(3);//Первая координата это время

//Начальные условия

X[0].push_back(0);

X[1].push_back(x0);

X[2].push_back(y0);

for (size_t k = 1; k < N; k++)

{

X[0].push_back(k*dt);

X[1].push_back

(

X[1].back() + dt * F(k*dt, X[1].back(), X[2].back())

);

X[2].push_back

(

X[2].back() + dt * G(k*dt, X[1].back(), X[2].back())

);

}

return X;

}

Метод Рунге для формулы Эйлера

std::vector> DSolve::EulerR(const int N, const double dt, std::function F, std::function G, double x0, double y0)

{

std::vector> XXN1 = Euler(N, dt, F, G, x0, y0);

std::vector> XXN2 = Euler(2*N, dt/2, F, G, x0, y0);

std::vector> XXN;

XXN.resize(3);

XXN[0].push_back(XXN1[0][0]);

XXN[1].push_back(

2 * XXN2[1][0] - XXN1[1][0]

);

XXN[2].push_back(

2 * XXN2[2][0] - XXN1[2][0]

);

for (size_t i = 1; i < N; i++)

{

XXN[0].push_back(XXN1[0][i]);

XXN[1].push_back(

2*XXN2[1][i*2]-XXN1[1][i]

);

XXN[2].push_back(

2 * XXN2[2][i*2] - XXN1[2][i]

);

}

return XXN;

}

Метод Адамса

std::vector> DSolve::Adams(const double N1, const double dt,

std::function F,

std::function G,

double x0, double y0)

{

const size_t N = N1;

std::vector> X;

X.resize(3);//Первая координата это время

//Начальные условия

X[0].push_back(0);

X[1].push_back(x0);

X[2].push_back(y0);

//N=1

X[0].push_back(1*dt);

std::vector> XX = RK4(2, dt, F, G, x0, y0);

X[1].push_back

(

XX[1].back()

);

X[2].push_back

(

XX[2].back()

);

for (size_t k = 1; k < N; k++)

{

X[0].push_back((k+1)*dt);

X[1].push_back

(

X[1][k] + (3.0/2)*dt * F(k*dt, X[1][k], X[2][k])

-

(1.0 / 2)*dt * F((k-1)*dt, X[1][k-1], X[2][k-1])

);

X[2].push_back

(

X[2][k] + (3.0/2)*dt * G(k*dt, X[1][k], X[2][k])

-

(1.0 / 2)*dt * G((k - 1)*dt, X[1][k - 1], X[2][k - 1])

);

}

return X;

}

Метод Рунге-Кутта 4 порядка

std::vector> DSolve::RK4(const double N1, const double dt,

std::function F,

std::function G,

double x0, double y0)

{

double k1, q1 = 0;

double k2, q2 = 0;

double k3, q3 = 0;

double k4, q4 = 0;

const size_t N = N1;

std::vector> X;

X.resize(3);//Первая координата это время

//Начальные условия

X[0].push_back(0);

X[1].push_back(x0);

X[2].push_back(y0);

for (size_t k = 1; k < N; k++)

{

k1 = F(k*dt, X[1].back(), X[2].back());

q1 = G(k*dt, X[1].back(), X[2].back());

k2 = F(k*dt+0.5*dt, X[1].back()+0.5*dt*k1, X[2].back()+0.5*dt*q1);

q2 = G(k*dt + 0.5*dt, X[1].back() + 0.5*dt*k1, X[2].back() + 0.5*dt*q1);

k3 = F(k*dt + 0.5*dt, X[1].back() + 0.5*dt*k2, X[2].back() + 0.5*dt*q2);

q3 = G(k*dt + 0.5*dt, X[1].back() + 0.5*dt*k2, X[2].back() + 0.5*dt*q2);

k4 = F(k*dt + dt, X[1].back() + dt*k3, X[2].back() + dt*q3);

q4 = G(k*dt + dt, X[1].back() + dt*k3, X[2].back() + dt*q3);

X[0].push_back(k*dt);

X[1].push_back

(

X[1].back() + (dt / 6.)*(k1 + 2 * k2 + 2 * k3 + k4)

);

X[2].push_back

(

X[2].back() + (dt / 6.)*(q1 + 2 * q2 + 2 * q3 + q4)

);

}

return X;

}

Конечно-разностный метод для решения краевой задачи в нормальной форме

std::vector BPSolve::FDM(const double h, std::function p, std::function q, std::function f, double a, double b, double A, double B, double a0, double a1, double b0, double b1)

{

std::vector X;

std::vector AA;

std::vector BB;

std::vector CC;

std::vector DD;

size_t N = b / h;

//x=0

AA.push_back(0);

BB.push_back(

a0*(1 - (h*p(0)) / 2)

+

a1 * ((q(0)*h) / 2 - 1 / h)

);

CC.push_back(

a1 / h

);

DD.push_back(

A*(1 - (h*p(0)) / 2)

+

a1 * f(0)*h / 2

);

for (size_t j = 1; j < N; j++)

{

AA.push_back(

1+p(j*h)/(2)

);

BB.push_back(

h*h*q(h*j)-2

);

CC.push_back(

1-(h*p(h*j)/(2))

);

DD.push_back(f(h*j)*h*h);

}

//x=b

BB.push_back(

b0*(1 + h * p(h*N) / 2)

+

b1 * (1 / h - (h*q(h*N)) / 2)

);

AA.push_back(

-b1 / h);

CC.push_back(0);

DD.push_back(

B*(1 + (h * p(h*N)) / 2)

-

(b1*h*f(h*N)) / 2

);

X = progonka(N, AA, BB, CC, DD);

return X;

}

Метод стрельбы для краевой задачи системы ОДУ

std::vector> BPSolve::SM(const double N1,

const double dt,

std::function F, std::function G,

double a, double b,

double A, double B,

double a0, double a1,

double b0, double b1)

{

DSolve solve;

double p1 = 0;

double p2 = 0;

double buff = 0;

double l1 = a;

double l2 = b;

const size_t N = N1;

//Выстрел

std::vector> Xl1;

std::vector> Xl2;

//Метод секущих

double l3 = 0;

do

{

Xl1 = solve.RK4(N1, dt, F, G, l1, A - (a0 / a1)*l1);

Xl2 = solve.RK4(N1, dt, F, G, l2, A - (a0 / a1)*l2);

p1 = b0 * Xl1[1].back() + b1 * Xl1[2].back()-B;

p2 = b0 * Xl2[1].back() + b1 * Xl2[2].back() -B;

l3 = l2 - ((l2 - l1)*p2) / (p2 - p1);

buff = abs(l3 - l2);

l1 = l2;

l2 = l3;
} while (buff > dt);

//Решаем краевую задачу

std::vector> X;

X = solve.RK4(N1, dt, F, G, l3, A - (a0 / a1)*l3);

return X;
}

1 2 3 4 5 6