Решение задач численными методами

ПОЯСНИТЕЛЬНАЯ ЗАПИСКА К курсовой работе по дисциплине «Информатика»

Тема Решение задач численными методами

Обычным дифференциальным уравнением n-го порядка называется выражение:

(1)

устанавливающее связь между независимой переменной x, искомой функцией y (x) и её производными до nго порядка включительно.

Решением дифференциального уравнения (1) называется функция y (x), которая при подстановке в уравнение (1) обращает его в тождество.

Так как таких решений у дифференциального уравнения в общем случае существует бесконечное множество, то каждое такое решение в отдельности называется частным решением дифференциального уравнения. Но все эти решения могут быть объединены в так называемое общее решение:

где — произвольные постоянные, принимающие для каждого отдельного частного решения отдельные значения, удовлетворяющие наперёд заданным значениям.

Если таким условием являются численно-заданные значения функции y (x) и её производные при одном и том же заданном значении аргумента, т. е.:

(2)

то совокупность дифференциального уравнения (1) с начальными условиями (2) называется задачей Коши интегрирования дифференциального уравнения (1).

Уравнение типа:

(3)

называется дифференциальным уравнением первого порядка.

Начальным условием для такого уравнения является задание значения функции в некотором значении аргумента:

(4)

Любое дифференциальное уравнение n-го порядка (1) может быть представлено системой n дифференциальных уравнений первого порядка:

(5)

Для решения последней используются методы, которые являются простым обобщением методов решения одного дифференциального уравнения первого порядка.

Так как подавляющее большинство дифференциальных уравнений решить аналитически очень трудно или невозможно, то огромное значение имеют методы их численного приближённого решения.

В данной курсовой производится численное решение задачи Коши для обыкновенного дифференциального уравнения первого и второго порядка методом Эйлера и Рунге — Кутты, и численное решение краевой задачи для ОДУ второго порядка.

Для расчетов используются пакет Mathsoft MathCAD и электронная таблица Microsoft Excel. Так же оба метода реализуются в программе, написанной на языке Microsoft Visual Basic.

дифференциальный эйлер mathcad

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

1.1 Аналитическое решение Для неоднородного уравнения первого порядка воспользуемся готовой формулой для задачи Коши:

Решение:

1.2 Решение в Mathcad

Численное решение задачи Коши для ОДУ первого порядка.

Задание. Найти решение задачи Коши с помощью метода Рунге-Кутты, разбив промежуток на 10 и 20 равных частей.

Построить графики полученных численных и аналитического решений.

Правая часть дифференциального уравнения Начальное условие Численное решение задачи Коши для системы ОДУ 1-го порядка (и одного ОДУ 1-го порядка) методом Рунге-Кутты находится с помощью стандартной функции rkfixed (y0,a, b, N, D)

y0 — вектор начальных условий для системы ОДУ 1-го порядка;

a, b — левый и правый конец промежутка, на котором ищется решение;

N — число разбиений промежутка [a;b];

D — вектор правых частей системы ОДУ 1-го порядка.

Для уравнения 1-го порядка вектор начальных условий и вектор правых частей имеют одну компоненту.

Функция возвращает матрицу, столбцами которой являются значения

1) независимой переменной x;

2) 1-й искомой функции;

и т.д.

Для уравнения 1-го порядка — одна искомая функция y.

Численное решение:

численное решение аналитическое решение значения x:

решение y:

Рис. 1

1.3 Решение в Excel

Численное решение задачи Коши для ОДУ первого порядка.

Задание. Найти решение задачи Коши с помощью методов Эйлера и Рунге-Кутты, разбив промежуток на 10 равных частей.

Построить графики полученных численных и аналитического решений.

a = 1,414 214

b = 2,414 214

y (2⁰, 5) = 1

Метод Эйлера:

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

i = 0,1,2,…, n-1 n= 10

где x[i] = x[0] + i*h

x[0] = a x[n] = b h=(b-a)/n = 0,1

Таблица 1

Метод Рунге-Кутты:

y[i+1] = y[i] + h*(k1+2*k2+2*k3+k4)/6

i = 0,1,2,…, n-1 n= 10

где k1 = f (x[i], y[i])

k2 = f (x[i]+h/2, y[i]+h*k½)

k3 = f (x[i]+h/2, y[i]+h*k2/2)

k4 = f (x[i]+h, y[i]+h*k3)

x[i] = x[0] + i*h

x[0] = a x[n] = b h=(b-a)/n = 0,1

Таблица 2

Рис. 2

1.4 Решение в Visual Basic

Окно программы:

Исходный код программы:

Private Function fun (X, Y)

y0 = -Y * X / (X ^ 2 + 1) + X

fun = y0

End Function

Private Sub Command1_Click ()

pr = 0

a = 2 ^ 0.5

b = 2 ^ 0.5 + 1

Y = 1

y00 = Y

n = Val (Text1.Text)

If n = 0 Then n = 10

If Option1. Value = True Then GoTo 1

If Option2. Value = True Then GoTo 2

If Option3. Value = True Then GoTo 3

Print #1, «Вычисление по схеме Эйлера»

Print #1,

Print #1, «i», «Xi», «Y1i», «f1()»

Print #1,

Print #1,

i = 0

X = a

f1 = fun (X, Y)

Y = Format (Y, «0.0000»)

X = Format (X, «0.0000»)

f1 = Format (f1, «0.0000»)

Print #1, i, X, Y, f1, f2

h = (b — a) / n

For X = (a + h) To b Step h

i = i + 1

Y = Y + h * f1

f1 = fun (X, Y)

Y = Format (Y, «0.0000»)

X = Format (X, «0.0000»)

f1 = Format (f1, «0.0000»)

Print #1, i, X, Y, f1

Next X

Print #1,

Print #1,

GoTo 4

h = (b — a) / n

Print #1, «Вычисление по схеме Рунге-Кутта»

Print #1,

Print #1, «i», «Xi», «Y1i», «K1», «K2», «K3», «K4»

Print #1,

i = 0

X = a

Y = y00

k1 = fun (X, Y)

k2x = X + h / 2

k2y = Y + (h / 2) * k1

k2 = fun (k2x, k2y)

k3x = X + h / 2

k3y = Y + (h / 2) * k2

k3 = fun (k3x, k3y)

k4x = X + h

k4y = Y + h * k3

k4 = fun (k4x, k4y)

Y = Format (Y, «0.0000»)

X = Format (X, «0.0000»)

k1 = Format (k1, «0.0000»)

k2 = Format (k2, «0.0000»)

k3 = Format (k3, «0.0000»)

k4 = Format (k4, «0.0000»)

k1 = fun (X, Y)

Print #1, i, X, Y, k1, k2, k3, k4

Y = Y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)

For X = (a + h) To b Step h

k1 = fun (X, Y)

k2x = X + h / 2

k2y = Y + (h / 2) * k1

k2 = fun (k2x, k2y)

k3x = X + h / 2

k3y = Y + (h / 2) * k2

k3 = fun (k3x, k3y)

k4x = X + h

k4y = Y + h * k3

k4 = fun (k4x, k4y)

Y = Format (Y, «0.0000»)

X = Format (X, «0.0000»)

k1 = Format (k1, «0.0000»)

k2 = Format (k2, «0.0000»)

k3 = Format (k3, «0.0000»)

k4 = Format (k4, «0.0000»)

i = i + 1

Print #1, i, X, Y, k1, k2, k3, k4

Y = Y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)

Next X

GoTo 4

Print #1, «Аналитический метод»

Print #1,

Print #1, «i», «Xi», «Yi»

Print #1,

Print #1,

i = 0

X = a

h = (b — a) / n

For X = a To b — h Step h

Y = (X ^ 2) / 3 + 1 / 3

Y = Format (Y, «0.0000»)

X = Format (X, «0.0000»)

Print #1, i, X, Y

i = i + 1

Next X

X = b

Y = (X ^ 2) / 3 + 1 / 3

Y = Format (Y, «0.0000»)

X = Format (X, «0.0000»)

Print #1, i + 1, X, Y

Print #1,

End Sub

Private Sub Form_Load ()

Open «c:coshi1.out» For Output As #1

End Sub

Private Sub Form_Unload (Cancel As Integer)

Close #1

End Sub

Распечатка файла Coshi1. out:

Таблица 3 Вычисление по схеме Эйлера Таблица 4 Вычисление по схеме Рунге-Кутта

Таблица 5 Аналитический метод дифференциальный эйлер mathcad

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

2.1 Решение в MathCad

Найти решение задачи Коши Правая часть дифференциального уравнения Начальное условие После введения новых неизвестных функций = = y' исходная задача сводится к задачи Коши для системы двух обыкновенных дифференциальных уравнений первого порядка Численное решение задачи Коши для системы ОДУ 1-го порядка (и одного ОДУ 1-го порядка) методом Рунге-Кутты находится с помощью стандартной функции rkfixed (y0,a, b, N, D)

y0 — вектор начальных условий для системы ОДУ 1-го порядка;

a, b — левый и правый конец промежутка, на котором ищется решение;

N — число разбиений промежутка [a;b];

D — вектор правых частей системы ОДУ 1-го порядка.

Функция возвращает матрицу, столбцами которой являются значения

1) независимой переменной x;

2) 1-й искомой функции;

3) 2-й искомой функции;

и т.д.

Численное решение:

y1(x) и y2(x) — значение первой и второй искомых функций значения x:

решение y:

Рис. 6

2.2 Решение в Excel

Задание. Найти приближенное решение задачи Коши.

a = 0,4

b = 1,4

y` = 1,4

y = 1,2

С помощью схем Эйлера и Рунге-Кутты при разбиении исходного промежутка на 10 равных частей.

Решение. Путем введения новых не известных функций y1 = y, y2=y0 задача Коши для обыкновенного дифференциального уравнения второго порядка сводится к задаче Коши для системы двух дифференциальных уравнений первого порядка Схема Эйлера применительно к задаче Коши для системы двух дифференциальных уравнений первого порядка

имеет вид:

Вычисления по схеме Эйлера помещаем в таблицу Таблица 6

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

где

k11 = f (x[i], y1[i], y2[i]);

k12 = f (x[i]+h/2, y1[i]+h/2*k11, y2[i]+h/2*k21);

k13 = f (x[i]+h/2, y1[i]+h/2*k12, y2[i]+h/2*k22);

k14 = f (x[i]+h, y1[i]+h*k13, y2[i]+h*k23);

k21 = g (x[i], y1[i], y2[i]);

k22 = g (x[i]+h/2, y1[i]+h/2*k11, y2[i]+h/2*k21);

k23 = g (x[i]+h/2, y1[i]+h/2*k12, y2[i]+h/2*k22);

k24 = g (x[i]+h, y1[i]+h*k13, y2[i]+h*k23);

Вычисления по схеме Рунге-Кутты оформляем в виде таблицы Таблица 7

Рис. 6

2.3 Решение в Visual Basic

Окно программы:

Исходный код программы:

Private Function fun (x, y1, y2)

y = 0.5 * x * y2 — y1 + 2

fun = y

End Function

Private Sub Command1_Click ()

pr = 0

a = 0.4

b = 1.4

y1 = 1.2

y2 = 1.4

n = Val (Text1.Text)

If n = 0 Then n = 10

If Option1. Value = True Then GoTo 1

If Option2. Value = True Then GoTo 2

Print #1, «Вычисление по схеме Эйлера»

Print #1,

Print #1, «i», «Xi», «Y1i», «Y2i», «f1()»; «f2()»

Print #1,

i = 0: x = a

f1 = y2: f2 = fun (x, y1, y2)

y1 = Format (y1, «0.0000»)

y2 = Format (y2, «0.0000»)

f1 = Format (f1, «0.0000»)

f2 = Format (f2, «0.0000»)

Print #1, i, x, y1, y2, f1, f2

h = (b — a) / n

For x = (a + h) To b Step h

i = i + 1

y1 = y1 + h * f1

y2 = y2 + h * f2

f1 = y2

f2 = fun (x, y1, y2)

y1 = Format (y1, «0.0000»)

y2 = Format (y2, «0.0000»)

f1 = Format (f1, «0.0000»)

f2 = Format (f2, «0.0000»)

Print #1, i, x, y1, y2, f1, f2

Next x

y1 = y1 + h * f1

y2 = y2 + h * f2

y1 = Format (y1, «0.0000»)

y2 = Format (y2, «0.0000»)

f1 = Format (f1, «0.0000»)

f2 = Format (f2, «0.0000»)

Print #1, i + 1, x, y1, y2

Print #1,

Print #1,

GoTo 3

h = (b — a) / n

Print #1, «Вычисление по схеме Рунге-Кутта»

Print #1,

Print #1, «i», «Xi», «Y1i», «Y2i», «K11», «K12», «K13», «K14», «K21», «K22», «K23», «K24»

Print #1,

i = 0

a = 0.4

b = 1.4

y1 = 1.2

y2 = 1.4

x = a

k11 = y2

k21 = fun (x, y1, y2)

k12 = y2 + h / 2 * k21

k22x = x + h / 2

k22y1 = y1 + (h / 2) * k11

k22y2 = y2 + (h / 2) * k21

k22 = fun (k22x, k22y1, k22y2)

k13 = y2 + h / 2 * k22

k23x = x + h / 2

k23y1 = y1 + (h / 2) * k12

k23y2 = y2 + (h / 2) * k22

k23 = fun (k23x, k23y1, k23y2)

k14 = y2 + h * k23

k24x = x + h

k24y1 = y1 + h * k13

k24y2 = y2 + h * k23

k24 = fun (k24x, k24y1, k24y2)

y1 = Format (y1, «0.0000»)

y2 = Format (y2, «0.0000»)

k11 = Format (k11, «0.0000»)

k12 = Format (k12, «0.0000»)

k13 = Format (k13, «0.0000»)

k14 = Format (k14, «0.0000»)

k21 = Format (k21, «0.0000»)

k22 = Format (k22, «0.0000»)

k23 = Format (k23, «0.0000»)

k24 = Format (k24, «0.0000»)

Print #1, i, x, y1, y2, k11, k12, k13, k14, k21, k22, k23, k24

y1 = y1 + h / 6 * (k11 + 2 * k12 + 2 * k13 + k14)

y2 = y2 + h / 6 * (k21 + 2 * k22 + 2 * k23 + k24)

k11 = y2

h = (b — a) / n

For x = (a + h) To b Step h

i = i + 1

k11 = y2

k21 = fun (x, y1, y2)

k12 = y2 + h / 2 * k21

k22x = x + h / 2

k22y1 = y1 + (h / 2) * k11

k22y2 = y2 + (h / 2) * k21

k22 = fun (k22x, k22y1, k22y2)

k13 = y2 + h / 2 * k22

k23x = x + h / 2

k23y1 = y1 + (h / 2) * k12

k23y2 = y2 + (h / 2) * k22

k23 = fun (k23x, k23y1, k23y2)

k14 = y2 + h * k23

k24x = x + h

k24y1 = y1 + h * k13

k24y2 = y2 + h * k23

k24 = fun (k24x, k24y1, k24y2)

y1 = Format (y1, «0.0000»)

y2 = Format (y2, «0.0000»)

k11 = Format (k11, «0.0000»)

k12 = Format (k12, «0.0000»)

k13 = Format (k13, «0.0000»)

k14 = Format (k14, «0.0000»)

k21 = Format (k21, «0.0000»)

k22 = Format (k22, «0.0000»)

k23 = Format (k23, «0.0000»)

k24 = Format (k24, «0.0000»)

Print #1, i, x, y1, y2, k11, k12, k13, k14, k21, k22, k23, k24

y1 = y1 + h / 6 * (k11 + 2 * k12 + 2 * k13 + k14)

y2 = y2 + h / 6 * (k21 + 2 * k22 + 2 * k23 + k24)

Next x

y1 = Format (y1, «0.0000»)

y2 = Format (y2, «0.0000»)

Print #1, i, x, y1, y2

End Sub

Private Sub Form_Load ()

Open «c:coshi2.out» For Output As #1

End Sub

Private Sub Form_Unload (Cancel As Integer)

Close #1

End Sub

Распечатка файла Coshi1. out:

Вычисление по схеме Эйлера Таблица 8

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

Найти численное решение краевой задачи

с помощью метода конечных разностей при разбиении исходного промежутка на 10, 20 и 40 равных частей.

Разбиение на 10 частей

Разбиение на 20 частей

Разбиение на 40 частей

Рис. 8

Заключение

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

1. Кудрявцев Е. М., «ДМК», Mathcad 8. Символьное и численное решение задач., «ДМК», 2010 г., 320 с.

2. Шафрин Ю. А. Основы компьютерной технологии. — М.: АБВ, 2009. — 656с.

3. Дьяконов В. Mathcad 2001: специальный справочник, СПб, 2011 г., 459 с.