Решение задач численными методами
ПОЯСНИТЕЛЬНАЯ ЗАПИСКА
К курсовой работе по дисциплине
"Информатика"
Тема Решение задач численными
методами
Введение
Обычным дифференциальным уравнением 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.2 Решение в Mathcad
Численное решение задачи Коши для ОДУ первого порядка.
Задание. Найти решение задачи Коши
с помощью метода Рунге-Кутты, разбив промежуток на 10 и 20 равных частей.
Построить графики полученных численных и аналитического решений.
Правая часть дифференциального уравнения
Начальное условие
Численное решение задачи Коши для системы ОДУ 1-го порядка (и одного ОДУ
1-го порядка) методом Рунге-Кутты находится с помощью стандартной функции rkfixed (y0,a,b,N,D)- вектор
начальных условий для системы ОДУ 1-го порядка;, b - левый и правый конец
промежутка, на котором ищется решение;- число разбиений промежутка [a;b];-
вектор правых частей системы ОДУ 1-го порядка.
Для уравнения 1-го порядка вектор начальных условий и вектор правых
частей имеют одну компоненту.
Функция возвращает матрицу, столбцами которой являются значения
) независимой переменной x;
) 1-й искомой функции;
и т.д.
Для уравнения 1-го порядка - одна искомая функция y.
Численное решение:
численное решение
аналитическое решение
значения x:
решение y:
Рис. 1
.3 Решение в Excel
Численное решение задачи Коши для ОДУ первого порядка.
Задание. Найти решение задачи Коши с помощью методов Эйлера и
Рунге-Кутты, разбив промежуток на 10 равных частей.
Построить графики полученных численных и аналитического решений.
a = 1,414214= 2,414214(2^0,5) = 1
Метод Эйлера:
[i+1] = y[i] + h*f(x[i],y[i])= 0,1,2,…,n-1 n= 10
где
x[i] = x[0] + i*h[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] ) = f( x[i]+h/2, y[i]+h*k1/2 )= f( x[i]+h/2, y[i]+h*k2/2 ) = f(
x[i]+h, y[i]+h*k3 )[i] = x[0] + i*h[0] = a x[n] = b h=(b-a)/n = 0,1
Таблица
2
Рис. 2
1.4 Решение в Visual Basic
Окно программы:
Исходный код программы:Function fun(X, Y)
= -Y * X / (X ^ 2 + 1) + X
= y0FunctionSub Command1_Click()= 0
a = 2 ^ 0.5= 2 ^ 0.5 + 1= 1= Y= Val(Text1.Text)n = 0 Then n =
10Option1.Value = True Then GoTo 1Option2.Value = True Then GoTo 2Option3.Value
= True Then GoTo 3 #1,
"Вычисление по схеме Эйлера"
Print #1,#1, " i", " Xi", "
Y1i", " f1()"#1,#1,= 0= a= fun(X, Y)= Format(Y,
"0.0000")= Format(X, "0.0000")= Format(f1,
"0.0000")#1, i, X, Y, f1, f2= (b - a) / nX = (a + h) To b Step h= i +
1= Y + h * f1= fun(X, Y)= Format(Y, "0.0000")= Format(X,
"0.0000")= Format(f1, "0.0000")#1, i, X, Y, f1X#1,#1,4= (b
- a) / n #1,
"Вычисление по схеме Рунге-Кутта"
Print #1,#1, " i", " Xi", "
Y1i", " K1", " K2", " K3", "
K4"#1,= 0= a= y00= fun(X, Y)x = X + h / 2y = Y + (h / 2) * k1= fun(k2x,
k2y)x = X + h / 2y = Y + (h / 2) * k2= fun(k3x, k3y)x = X + hy = Y + h * k3=
fun(k4x, k4y)= Format(Y, "0.0000")= Format(X, "0.0000")=
Format(k1, "0.0000")= Format(k2, "0.0000")= Format(k3,
"0.0000")= Format(k4, "0.0000")= fun(X, Y)#1, i, X, Y, k1,
k2, k3, k4= Y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)X = (a + h) To b Step h=
fun(X, Y)x = X + h / 2y = Y + (h / 2) * k1= fun(k2x, k2y)x = X + h / 2y = Y +
(h / 2) * k2= fun(k3x, k3y)x = X + hy = Y + h * k3= fun(k4x, k4y)= Format(Y,
"0.0000")= Format(X, "0.0000")= Format(k1,
"0.0000")= Format(k2, "0.0000")= Format(k3,
"0.0000")= Format(k4, "0.0000")= i + 1#1, i, X, Y, k1, k2,
k3, k4= Y + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)X 4
Print
#1, "Аналитический метод"
Print
#1,
Print #1, " i", " Xi", "
Yi"#1,#1,= 0= a= (b - a) / nX = a To b - h Step h= (X ^ 2) / 3 + 1 / 3=
Format(Y, "0.0000")= Format(X, "0.0000")#1, i, X, Y= i +
1X= b= (X ^ 2) / 3 + 1 / 3= Format(Y, "0.0000")= Format(X,
"0.0000")#1, i + 1, X, Y#1,SubSub Form_Load()"c:\coshi1.out"
For Output As #1SubSub Form_Unload(Cancel As Integer)#1Sub
Распечатка файла Coshi1.out:
Таблица 3 Вычисление по схеме Эйлера
Таблица 4 Вычисление по схеме Рунге-Кутта
Таблица 5 Аналитический метод
дифференциальный эйлер mathcad
2. Решение задачи Коши для ОДУ второго порядка
.1 Решение в MathCad
Найти решение задачи Коши
Правая часть дифференциального уравнения
Начальное условие
После введения новых неизвестных функций 
= 
= y' исходная задача сводится к
задачи Коши для системы двух обыкновенных дифференциальных уравнений первого
порядка
Численное решение задачи Коши для системы ОДУ 1-го порядка (и одного ОДУ
1-го порядка) методом Рунге-Кутты находится с помощью стандартной функции rkfixed(y0,a,b,N,D)- вектор
начальных условий для системы ОДУ 1-го порядка;, b - левый и правый конец
промежутка, на котором ищется решение;- число разбиений промежутка [a;b];-
вектор правых частей системы ОДУ 1-го порядка.
Функция возвращает матрицу, столбцами которой являются значения
) независимой переменной x;
) 1-й искомой функции;
и т.д.
Численное решение:
(x) и y2(x) - значение первой и второй искомых функций
значения x:
решение y:
Рис. 6
2.2 Решение в Excel
Задание. Найти приближенное решение задачи Коши.
= 0,4= 1,4` = 1,4= 1,2
С помощью схем Эйлера и Рунге-Кутты при разбиении исходного промежутка на
10 равных частей.
Решение. Путем введения новых не известных функций y1 = y, y2=y0 задача
Коши для обыкновенного дифференциального уравнения второго порядка сводится к
задаче Коши для системы двух дифференциальных уравнений первого порядка
Схема Эйлера применительно к задаче Коши для системы двух
дифференциальных уравнений первого порядка
имеет вид:
Вычисления по схеме Эйлера помещаем в таблицу
Таблица 6
Схема Рунге-Кутты применительно к задаче Коши для системы двух
дифференциальных уравнений первого порядка
где= f(
x[i], y1[i], y2[i] ); = f( x[i]+h/2, y1[i]+h/2*k11, y2[i]+h/2*k21 ); = f(
x[i]+h/2, y1[i]+h/2*k12, y2[i]+h/2*k22 ); = f( x[i]+h, y1[i]+h*k13, y2[i]+h*k23
); = g( x[i], y1[i], y2[i] );= g( x[i]+h/2, y1[i]+h/2*k11, y2[i]+h/2*k21 );= g(
x[i]+h/2, y1[i]+h/2*k12, y2[i]+h/2*k22 );= g( x[i]+h, y1[i]+h*k13, y2[i]+h*k23
);
Вычисления по схеме Рунге-Кутты оформляем в виде таблицы
Таблица 7
Рис. 6
2.3 Решение в Visual Basic
Окно программы:
Исходный код программы:
Function fun(x, y1, y2)= 0.5 * x * y2 - y1 + 2= yFunctionSub
Command1_Click()= 0= 0.4= 1.4= 1.2= 1.4= Val(Text1.Text)n = 0 Then n =
10Option1.Value = True Then GoTo 1Option2.Value = True Then GoTo 2 #1, "Вычисление по схеме
Эйлера"
Print #1,#1, " i", " Xi", "
Y1i", " Y2i", " f1()"; " f2()"#1,= 0: x = a=
y2: f2 = fun(x, y1, y2)= Format(y1, "0.0000")= Format(y2,
"0.0000")= Format(f1, "0.0000")= Format(f2,
"0.0000")#1, i, x, y1, y2, f1, f2= (b - a) / nx = (a + h) To b Step
h= i + 1= y1 + h * f1= y2 + h * f2= y2= fun(x, y1, y2)= Format(y1,
"0.0000")= Format(y2, "0.0000")= Format(f1,
"0.0000")= Format(f2, "0.0000")#1, i, x, y1, y2, f1, f2x=
y1 + h * f1= y2 + h * f2= Format(y1, "0.0000")= Format(y2,
"0.0000")= Format(f1, "0.0000")= Format(f2,
"0.0000")#1, i + 1, x, y1, y2#1,#1,3= (b - a) / n #1, "Вычисление по схеме
Рунге-Кутта"
Print #1,#1, " i", " Xi", "
Y1i", " Y2i", " K11", " K12", "
K13", " K14", " K21", " K22", "
K23", " K24"#1,= 0= 0.4= 1.4= 1.2= 1.4= a= y2= fun(x, y1, y2)=
y2 + h / 2 * k21x = x + h / 2y1 = y1 + (h / 2) * k11y2 = y2 + (h / 2) * k21=
fun(k22x, k22y1, k22y2)= y2 + h / 2 * k22x = x + h / 2y1 = y1 + (h / 2) * k12y2
= y2 + (h / 2) * k22= fun(k23x, k23y1, k23y2)= y2 + h * k23x = x + hy1 = y1 + h
* k13y2 = y2 + h * k23= fun(k24x, k24y1, k24y2)= Format(y1,
"0.0000")= Format(y2, "0.0000")
k11 = Format(k11, "0.0000")= Format(k12,
"0.0000")= Format(k13, "0.0000")= Format(k14,
"0.0000")= Format(k21, "0.0000")= Format(k22,
"0.0000")= Format(k23, "0.0000")= Format(k24,
"0.0000")
Print #1, i, x, y1, y2, k11, k12, k13, k14, k21, k22, k23,
k24= y1 + h / 6 * (k11 + 2 * k12 + 2 * k13 + k14)= y2 + h / 6 * (k21 + 2 * k22
+ 2 * k23 + k24)= y2= (b - a) / nx = (a + h) To b Step h= i + 1= y2= fun(x, y1,
y2)= y2 + h / 2 * k21x = x + h / 2y1 = y1 + (h / 2) * k11y2 = y2 + (h / 2) *
k21= fun(k22x, k22y1, k22y2)= y2 + h / 2 * k22x = x + h / 2y1 = y1 + (h / 2) *
k12y2 = y2 + (h / 2) * k22= fun(k23x, k23y1, k23y2)= y2 + h * k23x = x + hy1 =
y1 + h * k13y2 = y2 + h * k23= fun(k24x, k24y1, k24y2)= Format(y1,
"0.0000")= Format(y2, "0.0000")
k11 = Format(k11, "0.0000")= Format(k12,
"0.0000")= Format(k13, "0.0000")= Format(k14,
"0.0000")= Format(k21, "0.0000")= Format(k22,
"0.0000")= Format(k23, "0.0000")= Format(k24,
"0.0000")
Print #1, i, x, y1, y2, k11, k12, k13, k14, k21, k22, k23,
k24= y1 + h / 6 * (k11 + 2 * k12 + 2 * k13 + k14)= y2 + h / 6 * (k21 + 2 * k22
+ 2 * k23 + k24)x= Format(y1, "0.0000")= Format(y2,
"0.0000")#1, i, x, y1, y2SubSub Form_Load()"c:\coshi2.out"
For Output As #1SubSub Form_Unload(Cancel As Integer)#1Sub
Распечатка файла Coshi1.out:
Вычисление по схеме Эйлера
Таблица 8
3. Численное решение краевой задачи для ОДУ второго порядка
Найти численное решение краевой задачи
с помощью метода конечных разностей при разбиении исходного промежутка на
10, 20 и 40 равных частей.
Разбиение на 10 частей
Разбиение на 20 частей
Разбиение на 40 частей
Рис. 8
Заключение
В проделанной курсовой работе все расчёты произведены с помощью численных
методов решения математических задач, которые были изучены ранее. В качестве
дополнительных материалов были использованы лабораторные работы и конспект
лекций по предмету.
Список литературы
1. Кудрявцев Е.М., "ДМК", Mathcad 8. Символьное и численное решение задач.,
"ДМК", 2010 г., 320 с.
. Шафрин Ю.А. Основы компьютерной технологии. - М.: АБВ,
2009. - 656с.
. Дьяконов В. Mathcad 2001: специальный справочник, СПб,
2011г., 459 с.