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

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

Курсовая робота

Методы приближеного решения дифференциальных уравнений

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

1.1 Постановка задания

Пусть задано обыкновенное дифференциальное уравнение

y(x) f (x, y(x)), x[a, b], (1)

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

где f(x) известная функция, y(x) искомая функция на [a, b], а также начальное условие:

y(x₀) y₀, x₀[a, b]. (2)

1. Разработать подпрограммы решения задачи Коши (1), (2) методом Эйлера, методом Рунге-Кутта порядка точности s и явным s-шаговым методом Адамса (s >1 указано в варианте задания). Разработанная программа должна давать численное решение задачи Коши каждым из упомянутых методов с шагом интегрирования

,

где N заданное натуральное число.

2. Найти решение задачи Коши (1), (2) в квадратурах.

3. Распечатать в виде сравнительной таблицы значения приближенных решений, полученных в соответствии с п. 1, а также вычисленные на ЭВМ значения точного решения y(x) в тех же точках.

4. Построить и сравнить между собой графики полученных решений.

Функция f(x), значения a, b, x₀, y₀ определяются вариантом задания.

Индивидуальное задание:

;

y(0)=0;

1.2 Основные теоретические сведения

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

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

Задача нахождения решения уравнения, удовлетворяющего заданным начальным условиям, называется задачей Коши.

Известно, что если функция f непрерывна по всем своим аргументам в некоторой замкнутой области и удовлетворяет условию Липшица по переменным , то задача Коши имеет единственное решение.

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

Итак, пусть на отрезке рассматривается следующая задача Коши:

(1)

(2)

Будем искать приближенные значения ее решения лишь в фиксированных точках Данного отрезка. Точки будем считать равноотстоящими:

по формуле Ньютона-Лейбница

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

(3)

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

Метод Эйлера

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



Отбросив член порядка и обозначив через приближенное значение , получим формулу

(4)

по которой, зная значение , можно последовательно вычислить приближенные значения решения во всех точках

Численный метод решения задачи (1), (2), определяемый формулой (4), называется методом Эйлера.

Метод Эйлера является простейшим численным методом интегрирования дифференциального уравнения. Недостатком этого метода является малая точность (из-за малой точности квадратурной формулы левых прямоугольников) и систематическое накопление ошибок.

Метод Эйлера относится к классу одношаговых методов. Это такие методы решения задачи (1), (2), которые позволяют найти приближенное значение решения этой задачи в узле по информации об этом решении в одном предыдущем узле .

Среди одношаговых методов можно выделить явные и неявные. В достаточно общем виде явные одношаговые методы можно записать в виде

Ясно, что метод Эйлера относится к классу явных методов.

Общий вид неявных одношаговых методов

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

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

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

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

Введем три набора параметров

выбором которых распорядимся в дальнейшем. С помощью параметров и составим величины



Если параметры и выбраны, то значения величин вычисляются последовательно.

Заметим, что величины

вообще говоря, не равны

Однако, при удачном выборе параметров можно считать, что такое равенство имеет место приближенно. Поэтому можно надеяться с помощью параметров (P) составить такую линейную комбинацию величин , которая будет являться аналогом квадратурных сумм и позволит найти приближенное значение интеграла в (3):

(5)

(здесь для упрощения записи опущен индекс n, означающий номер шага).

Рассмотрим теперь проблему выбора введенных параметров. Величина

представляет собой погрешность приближенного равенства (5). Если правая часть уравнения (1) - достаточно гладкая функция, то функция будет иметь достаточно большое число производных. Тогда по формуле Маклорена

(6)

Если удастся подобрать параметры так, чтобы при j=0,1,…, s, а , то погрешность формулы (5) будет величиной порядка :

(7)

Число S называют порядком (или степенью) точности данного метода типа Рунге-Кутта. Это название связано с тем, что если на каждом шаге погрешность метода имеет порядок , то погрешность на всем рассматриваемом отрезке будет величиной порядка .

Для построения по методу Рунге-Кутта при данном r одношаговых правил возможно более высокого порядка точности S выражают величины через значения функции и ее частных производных в точке (x, y) и требуют, чтобы для любой гладкой функции f (x, y) обращалось в нуль возможно большее число этих величин. Иными словами, параметры выбирают из требования, чтобы разложение

(8)

и разложение линейной комбинации по степеням h совпадали для произвольной функции f (x, y) для членов с возможно более высокими степенями h [2,3].

Записать в общем виде систему уравнений для определения параметров трудно. Поэтому рассмотрим лишь несколько примеров построения одношаговых методов по способу Рунге-Кутта.

Экстраполяционный метод Адамса

Экстраполяционный метод Адамса относится к явным численным методам решения задачи Коши.

В основе большинства численных методов лежит формула (3)

Будем вычислять этот интеграл путем замены подинтегральной функции ее интерполяционным многочленом, построеннном по предшествующим (к+1) точкам . Так как эти точки не принадлежат промежутку интегрирования , то фактически это будет экстраполирование. Отсюда название метода.

Поскольку узлы интерполирования равноотстоящие, а отрезок, на котором требуется аппроксимировать функцию , расположен справа от таблицы узлов, то в качестве интерполяционного многочлена естественно выбрать многочлен Ньютона для интерполирования назад. В формуле (3) сделаем замену переменных , в результате которой эта формула примет вид

(4)

и подставим функцию как сумму ее интерполяционного многочлена и остаточного члена



(5)

(6)

Подставим в (4) под знак интеграла вместо ее представление по формуле (5). Получим:

(7),

где ,

Формула (7) справедлива для всех значений .

Рассмотрим остаточный член формулы (7):

Поскольку сохраняет знак на [0,1], а функция - достаточно гладкая, то можно воспользоваться обобщенной теоремой о среднем значении.

Или, если обозначить , то

- ограниченная величина. (8)

Если h достаточно мало, то R_k мала и, отбрасывая R_k в формуле (7), и переходя к обозначению получим следующую расчетную формулу экстраполяционного метода Адамса:

, (9),

где .

Формула (9) является явной, при чем для вычисления используются - (k+1) предыдущих значений, то есть метод является (k+1) - шаговым.

Вычисления по формуле (9) начинаются с n=k, следовательно, необходимо до начала вычислений по формуле (9) знать предварительно значения (их можно посчитать любым другим одношаговым численным методом желательно того же порядка точности).

Минимальная погрешность формулы (9) на основании оценки (8) =, а глобольная погрешность =.

Замечание:

Всякий S-шаговый метод Адамса имеет на отрезке [x₀, X] наивысший порядок точности S.

Правило Рунге приближенной оценки погрешности

Рассмотрим на отрезке (х₀, Х) задачу Коши (1), (2).

(1)

у(х₀)=у₀(2)

Выполним разбиение этого отрезка на N равных частей с шагом h точками х_n=x₀+n*h, n=0,…, N-1 и будем применять для приближенного решения задачи (1), (2) некоторый явный одношаговый метод S-того порядка точности. Обозначим через у(х_n) - значение точного решения в т.х_n. y_n - приближенное значение решения в т. х_n ().

Напомним, что численный метод решения задачи Коши имеет S-тый порядок точности, если существует такое число S>0, что при достаточно малом h (y(x_n) - y_n)=O(h^s) для любых точек из отрезка [х₀, Х]). Из этой формулы в частности выплывает, что

y(x_n) - y_nC*h^s, (3)

где С является константой, зависящая вообще говоря от x_n.

Предположим, что в некоторой точке данным одношаговым методом S-того порядка точности получены приближенные решения задачи Коши с шагами h₁, h₂. Обозначим соответствующие приближенные решения через y_h1, y_h2.

Тогда на основании формулы (3) можно записать

(4),

(5),

Вычтем из (5) (4)

(6),

Подставляя значения C(t) в формулу (4) найдем приближенное значение для погрешности решения y_h1:

(7),

Подставляя в формулу (5) получим приближенное значение для y_h2:

(8),

Формула (7) носит название правило Рунге приближенной оценки погрешности (первое правило Рунге). В частности, если h₁=h, а h₂=2h, то тогда формула (7) примет вид:

Выразим из формулы (7) и (8) y(f) получим:

(9),

Формула (9) называется правилом Рунге уточнения результата (второе правило Рунге). Следует ожидать. что приближенное значение решения в правой части формулы(9) будет точнее, чем каждый из приближенных результатов y_h2 и y_h1 полученных в каждой точке t.

1.3 Численный эксперимент и анализ результатов

По заданию нужно решить конкретную задачу Коши:

;

у(0)=0

После запуска программы были получены следующие результаты:

Метод Рунге-Кутта		Метод Адамса		Метод Эйлера
X	Y		x	y		x	y
0	0		0	0		0	0
0,1	0,02		0,1	0,02		0,1	0,02
0,2	0,04		0,2	0,04		0,2	0,04
0,3	0,06		0,3	0,06		0,3	0,06
0,4	0,08		0,4	0,08		0,4	0,08
0,5	0,1		0,5	0,1		0,5	0,1
0,6	0,12		0,6	0,12		0,6	0,12
0,7	0,14		0,7	0,14		0,7	0,14
0,8	0,16		0,8	0,16		0,8	0,16
0,9	0,18		0,9	0,18		0,9	0,18
1	0,2		1	0,2		1	0,2
Метод Рунге-Кутта для 2h	Метод Адамса для 2h	Метод Эйлера для 2h
X	Y		x	y		x	y
0	0		0	0		0	0
0,2	0,04		0,2	0,04		0,2	0,04
0,4	0,08		0,4	0,08		0,4	0,08
0,6	0,12		0,6	0,12		0,6	0,12
0,8	0,16		0,8	0,16		0,8	0,16
1	0,2		1	0,2		1	0,2
погр. Рунге-Кутта		погр. Адамса		погр. Эйлера
0	0		0	0		0	0
0,2	0		0,2	0		0,2	0
0,4	0		0,4	0		0,4	0
0,6	0		0,6	0		0,6	0
0,8	0		0,8	0		0,8	0
1	0		1	0		1	0

График решения имеет вид:

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

Оценка по первому правилу Рунге показала нулевые погрешности.

2. Методы приближенного решения граничных задач для ОДУ

2.1 Постановка задания

Пусть задано линейное дифференциальное уравнение

y» (x) + p(x) y(x) + q(x) y(x) f(x), x[a, b] (1)

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

₀ y(a) + ₁ y(a) A; ₀ y(b) + ₁ y(b) B. (2)

Функции p(x), q(x), f(x) предполагаются известными и непрерывными на [a, b]. Числовые коэффициенты ₀, ₁, А, ₀, ₁, В также известны, причем для коэффициентов ₀, ₁, ₀, ₁ выполняются условия ₀+₁ 0, ₀+₁ 0.

Найти решение y(x), x[a, b] задачи (1), (2). С этой целью:

1. Разработать подпрограммы приближенного решения краевой задачи (1), (2) двумя методами: численным (методом конечных разностей) и аналитическим (методом коллокаций).

2. Найти решение конкретного варианта краевой задачи (1), (2) в квадратурах.

3. Распечатать в виде сравнительной таблицы значения приближенных решений, полученных указанными методами в точках отрезка [a, b]

, , , N>0,

а также вычисленные на ЭВМ значения точного решения y(x) в тех же точках.

4. Построить и сравнить между собой графики полученных решений.

5. Провести анализ результатов.

Замечание:

1) Если граничные условия (2) являются условиями первого рода:

y(a) A, y(b) B, (2)

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

, ; (3)

, .

2) Если не оговорено противное, то рекомендуется принять N=10, n=3.

Индивидуальное задание:

;

2.2 Основные теоретические сведения

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

Такие задачи называются многоточечными.

Математическая постановка задачи. Пусть на отрезке (a, b) задано обычное дифференциальное уравнение n-того порядка:

y⁽ⁿ⁾(x)=f (x, y^I(x), y^II(x), …, y^(n-1)(x)).

Пусть имеются k точек таких, что a и n соотношений связывающих значения искомого решения y(x) и производных до (n-1) - го порядка включительно в этих точках, а именно:

(2),

Здесь y, Q_j известные непрерывные функции своих аргументов. А_j - известные константы.

Задача нахождения решения уравнения (1) удовлетворяющего условиям (2) называется k-точечной задачей. При этом в частности, если k=1 и х₁ совпадает с точкой а, (х₁=b), то задача 1,2 называется задачей Коши. Если k=2 и х₁=а, х₁=b, то задача называется краевой (граничной задачей). В дальнейшем будем всегда предполагать, что решение задачи (1), (2) существует, единственно и обладает достаточной степенью гладкости.

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

Краевая задача называется линейной, если линейно дифференциальное уравнение и дополнительные условия.

Известно, что краевая задача могут быть решены точно лишь в крайне редких частных случаях, поэтому в основном они решаются приближенно. Все приближенные методы решения краевых задач делятся на 3 группы:

1) аналитические методы

а) метод колокации

б) метод Галёркина

в) обобщенный метод моментов и др.

2) численные методы

а) конечных разностей

б) конечных элементов и др.

3) методы численно-аналитические

а) методы редукции линейно краевой задачи к двум задачам Коши

б) метод стрельбы и др.

Метод конечных разностей (метод сеток)

Метод сеток или метод конечных разностей - это численный метод решения краевых задач, наиболее универсальный из всех приближенных методов (применим к линейным и нелинейным краевым задачам, применим к задачам с высокой размерностью, для областей со сложной формы).

Основная идея метода конечных разностей состоит в сведении краевой задачи к решению некоторой системы алгебраических уравнений.

Идею метода сеток рассмотрим на примере решения линейной краевой задачи для обыкновенного дифференциального уравнения 2 - го порядка.

Постановка задачи:

(1)

Граничное условие:

(2)

(3)

Где f, p(x), q(x) - известны, непрерывны функции на (a, b)

- известные константы.

у(х) - искомая функция на [a, b].

Разобьем отрезок (а, b) на N > 0 равных для простоты частей, точками

,

Совокупность указанных точек называется сеткой на отрезке [a, b], а каждая из этих точек называется узлом сетки , .

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

Поставим задачу выразить значения производных через значения самой функции в точках сетки, для этого запишем уравнение (1) во внутренних узлах сетки:

(4)

Запишем разложение функции в ряд Тейлора в окрестности точки :

Рассмотрим для точки x=x_k+1:

(5) у(x_k+1)=y(x_k) + y`(x_k)*h + y``(x_k) + y```(x_k) +O()

А также для точки x=x_k-1:

(6) у(x_k-1)=y(x_k) - y`(x_k)*h + y``(x_k) - y```(x_k) +O()

Если выразить из (5) и (6) значение у` в точке x_k, то получим следующие приближённые выражения:

(7) у`(x_k) ? - правая разностная производная функции у(х) в точке x_k,

Погрешность = О(h).

А также из (6) получим:

(7) у`(x_k) ? - левая разностная производная функции у(х) в точке x_k,

Погрешность = О(h).

Теперь вычтем из (5) формулу (6). Получим:

(9) у`(x_k) ? - центральная разностная производная функции у(х) в точке x_k,

Погрешность = О(h²).

Сложим формулы (5) и (6):

(10) у``(x_k) ? - разностная производная второго порядка функции у(х) в точке x_k,

Подставим в (4) вместо у` и у`` их выражения по формулам (9) и (10).Получим:

+ p(x_k)* + q(x_k)*y(x_k) ? f(x_k), k=1, n-1.

Умножим обе части на , приведём к обозначению у_k ? у(x_k):

(11) (1+*p(x_k)) *y_k+1 + (q(x_k) *h² - 2)*y_k + (1 - *p(x_k)) *y_k-1 = f(x_k)*h².

Соотношение (11) является разностным аналогом ДУ (1) и аппроксимирует его с погрешностью О(h²).

Построим теперь разностные аналоги граничных условий. Заменим:

? у`(а) во (2) правой разностной производной по формуле (7);

? у`(b) в (3) левой разностной производной по формуле (8);

Из второго граничного условия и формулы (7) будем иметь:

Из (2) и (7): б ₀ *y(a) + б ₁ * ? A.

Из (3) и (8): в ₀ *y(a) + в ₁ * ? B.

Домножим на h, приведём подобные слоагаемые, у(а) =у₀ вынесем за скобки:

(12) (б ₀h- б ₁) у₀ + б ₁у₁ = Ah.

(13) (в ₀h- в ₁) у_n + в ₁у_n-1 = Bh.

Формулы (12), (13) являются разностными аналогами граничных условий (2) и (3) соответственно, и аппроксимируют их с погрешностью О(h²).

Таким образом, соотношение (11), (12), (13) образуют разностную схему для исходной задачи (1), (2), (3), и представляют собой СЛАУ относительно приближённых значений у₀, у₁, …, у_n искомого решения у(x) в точках сетки. То есть наша задача свелась к решению СЛАУ с (N+1) неизвестной и (N+1) уравнением.

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

Кроме того, это система обладает трёхдиагональной матрицей:

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

Метод разностной прогонки

Метод разностной прогонки - это прямой метод решения СЛАУ с матрицей трехдиагонального вида. Он является разновидностью метода Гаусса исключения неизвестных.

Рассмотрим систему сеточных уравнений (11), (12), (13). Предположим, что в уравнениях (11) при любом k= .

Для удобства разделим на него каждое из уравнений (11) и перепишем их в следующем виде:

; k= (11')

Из уравнения (12) выразим y₀. Получим:

Введем обозначения:

, (14)

Тогда (15)

Основная идея метода разностной прогонки состоит в том, чтобы связать любые 2 соседних значения y_k, y_k+1 соотношениями вида (15), то есть чтобы получить формулу

, k=, (16)

где C_k и d_k - неизвестные пока коэффициенты.

То есть граничное условие (12) как бы прогоняется через все искомые значения y_k.

Присоединим формулу (15) к соотношению (16)

, k=, (17)

С помощью формулы (17) исключим y_k-1 из уравнения (11'):

, k=

, k=

Сравнивая полученные соотношения с формулами (16) находим прогонные коэффициенты:

; . (18)

Зная значения C₀ и d₀ по формулам(14) можно последовательно отыскать C_k и d_k по формулам (18).

Процесс вычисления прогоночных коэффициентов называется прямым ходом метода разностной прогонки.

Из формул (16) найдем при k=N-1:

Подставим это значение в формулу (13):

(19)

С помощью формулы (19) и (17) можно последовательно отыскать интересующие нас значения y_N, y_N-1, y_N-2,…, y₁, y₀.

Замечание:

Если в исходной системе (11), (12), (13) диагональные элементы доминируют по строкам или столбцам, то метод разностной прогонки обладает корректным и вычислительно устойчивым алгоритмом.

Метод коллокаций

Пусть на [a, b] задана ЛКЗ для ОДУ второго порядка:

(1)

Где p(x), q(x) f(x) - известные непрерывные на (а, b) функции.

Имеются линейные граничные условия:

(2)

(3)

Где f, p(x), q(x) - известные, непрерывные функции на (a, b)

- известные константы.

у(х) - искомая функция на [a, b].

Введем для краткости следующее обозначение:

где L - линейный оператор:

.

Тогда (1) можно переписать:

Ly(x)=f(x), x є (a, b)

Будем искать решение задачи (1), (2), (3) в следующем виде:

(5) , для всех х є [a, b]

где n є N - натуральное число, выбирается из соображений точности

- известные базисные функции, на которые накладываются ограничения.

k=1, n - числовые коэфициенты, подлежащие нахождению.

Требования к базисным функциям:

1. - линейно независимые на [a, b].

2.

- образуют полную систему функций в пространстве на отрезке [a, b]. То есть любую функцию y(x) из этого пространства можно сколь угодно точно приблизить линейной комбинацией этих базисных функций.

3. удовлетворяют соответствующим однородным граничным условиям:

удовлетворяет условиям:

(2)

(3)

4. - дважды непрерывно дифференцируемые на [a, b].

Из условия 3. вытекает что функция вида (5) будет удовлетворять граничным условиям (2) и (3) при любых значениях коэффициентов .

Однако функция не обязана удовлетворять ДУ (1) для любого х из [a, b]:

Введем в рассмотрение систему точек , которые назовем точками коллокаций и потребуем чтобы уравнение Ly(x)=f(x), x є (a, b) для функции выполнялось в этих точках за счет выбора коэффициентов и потребуем:

(6)

Из (6) получаем:

(7)

(7) - система линейных алгебраических уравнений относительно неизвестных а1,…, аn.

Матрица этой системы имеет следующий вид:



Решив систему (7) любым известным способом найдем коэффициенты и подставив в формулу (5) получим решение исходной краевой задачи.

Замечание

В теории метода доказано что det СЛАУ 0 при достаточно больших n.

Обусловленность СЛАУ. Число обусловленности матрицы

Пусть задана система линейных алгебраических уравнений , , где - это квадратная матрица n - го порядка.

Обозначим через обратную матрицу.

Определение:

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

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

Обусловленность матрицы характеризуется следующим числом

,

которое называется числом обусловленности матрицы в рассматриваемой норме.

Чем больше число обусловленности матрицы, тем хуже она обусловлена.

Если число обусловленности ближе к 1, то матрица хорошо обусловлена, если оно порядка , то матрица является плохо обусловленной.

2.3 Численный эксперимент и анализ результатов

;

Из заданного ДУ легко видеть что:

Метод коллокации

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

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

Вычислим значения точек коллокации:

x[1] = 0,0670

x[2] = 0,5000

x[3] = 0,9330

Подставив найденные значения точек коллокации в базисные функции, получим следующую систему линейных алгебраических уравнений:

-2,0580 1,5900 0,3473 = 0,9330

-2,0000 -1,1250 -0,1250 = 0,5000

-1,1920 -2,9025 -4,2535 = 0,0670

Решив систему относительно неизвестных , получим:

a[1] = -0,3367

a[2] = 0,1573

a[3] = -0,0287

Вычислим число обусловленности данной системы уравнений. Для этого запишем систему в матричном виде Ах=f

-2,0580 1,5900 0,3473 a₁ 0,9330

А= -2,0000 -1,1250 -0,1250; x= a_2; f = 0,5000

-1,1920 -2,9025 -4,2535 a₃ 0,0670

Число обусловленности:

C(A) = ||А||₂ *||(А)^-1||₂

||А||₂ = 6,7258

	-0,2122	-0,27615	-0,00921
(А)-1=	0,401053	-0,43991	0,045674
	-0,2142	0,377571	-0,26369

||(А)^-1||₂=1,093631

C(A) = ||А||2 *||(А) - 1||2= 6,7258*1,093631= 6,956586791

Можно говорить, что матрица хорошо обусловлена.

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

Вычисление значений функции на отрезке [0; 0,1] с шагом 0,1

Y (0,000) = 2,0000

Y (0,100) = 2,0711

Y (0,200) = 2,1510

Y (0,300) = 2,2387

Y (0,400) = 2,3332

Y (0,500) = 2,4337

Y (0,600) = 2,5394

Y (0,700) = 2,6495

Y (0,800) = 2,7633

Y (0,900) = 2,8803

Y (1,000) = 3,0000

Сводная таблица значений приближенного решения краевой задачи

X	коллокации	конечных разностей
0	2	2
0,1	2,0711	2,0668
0,2	2,151	2,1436
0,3	2,2387	2,2292
0,4	2,3332	2,3224
0,5	2,4337	2,4225
0,6	2,5394	2,5286
0,7	2,6495	2,64
0,8	2,7633	2,756
0,9	2,8803	2,8762
1	3	3



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

Перечень использованной литературы

1. Балашова С.Д. Чисельн методи: Ч. 2. Методи розв'язування диференцальних та нтегральних рвнянь: Навч. посбник. - К.: НМК ВО, 1992.

2. Борис Павлович Демадозич и Исаак Абрамович Марон «Основы вычислительной математики», М., 1966 р. 664 стр., Издательство «Наука».

Приложения

Приложение 1. Программа решения задачи Коши

#include <stdio.h>

#include <iostream>

#include <stdlib.h>

#include <conio.h>

#include <math.h>

#include <malloc.h>

#include <windows.h>

using namespace std;

double s4itat_f (double xi, double yi)

{

double f;

f=(0.5)*xi*yi - (2.5)*yi*yi+(0.2);

return f;

}

void nayti_to4ki (double *x, int N, double h)

{for (int n=1; n<N; n++)

{

x[n]=x[0]+n*h;

}

}

void rewit_eyler (double *x, double *y, int N, double h)

{

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

{

y [n+1]=y[n]+h*s4itat_f (x[n], y[n]);

}

}

void rewit_rungekutt (double *x, double *y, int N, double h, int r)

{

double **K;

K=(double **) malloc (sizeof(double *)*(N+1));

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

{

K[n]=(double *) malloc (sizeof(double)*r);

}

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

{

K[n] [1-1]=h*s4itat_f (x[n], y[n]);

K[n] [2-1]=h*s4itat_f((x[n]+h/2), (y[n]+(K[n] [1-1])/2));

K[n] [3-1]=h*s4itat_f((x[n]+h/2), (y[n]+(K[n] [2-1])/2));

K[n] [4-1]=h*s4itat_f (x[n+1], (y[n]+K[n] [3-1]));

y [n+1]=y[n]+(K[n] [1-1]+2*K[n] [2-1]+2*K[n] [3-1]+K[n] [4-1])/6;

}

free(K);

}

void rewit_adams (double *x, double *y, int N, double h, int r)

{

for (int n=r-1; n<N; n++)

{

//y [n+1]=y[n]+h*(23*s4itat_f (x[n], y[n]) - 16*s4itat_f (x[n-1], y [n-1])+5*s4itat_f (x[n-2], y [n-2]))/12;

y [n+1]=y[n]+h*(55*s4itat_f (x[n], y[n]) - 59*s4itat_f (x[n-1], y [n-1])+37*s4itat_f (x[n-2], y [n-2]) - 9*s4itat_f (x[n-3], y [n-3]))/24;

}

}

void vivesti_rezultati (double *x, double *y_eyler, double *y_rungekutt, double *y_adams, int N, FILE *f)

{

setlocale (0, "»);

char s[1250];

printf (» -\n»);

printf (» |Точки |Эйлера |Рунге-Кутта |Адамса |\n»);

printf (» -\n»);

fprintf (f, «|Точки |Эйлера |Рунге-Кутта |Адамса |\n»);

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

{

printf (» |%f\t|%f\t|%f\t|%f\t|\n», x[n], y_eyler[n], y_rungekutt[n], y_adams[n]);

fprintf (f, «|%f\t|%f\t|%f\t|%f\t|\n», x[n], y_eyler[n], y_rungekutt[n], y_adams[n]);

}

cout<<» -\n»;

cout<<endl;

}

int main()

{

double *y_eyler, *y_rungekutt, *y_adams, *x;

double h;

int N;

int r;

FILE *f;

f=fopen («out.txt», «w»);

setlocale (0, "»);

cout<< «Количество точек N=»;

cin>>N;

//N=10;

r=4;

y_eyler=(double *) malloc (sizeof(double)*(N+1));

y_rungekutt=(double *) malloc (sizeof(double)*(N+1));

y_adams=(double *) malloc (sizeof(double)*(N+1));

x=(double *) malloc (sizeof(double)*(N+1));

y_eyler[0]=0; // начальное условие

y_rungekutt[0]=0;

y_adams[0]=0;

x[0]=0;

x[N]=1;

h=(x[N] - x[0])/N;

//h=0.1;

cout<< «f=0.5*xy-2.5*y^2+0.2»<<endl;

cout<< «N=»<<N<<endl;

cout<< «h=»<<h<<endl;

nayti_to4ki (x, N, h);

rewit_eyler (x, y_eyler, N, h);

rewit_rungekutt (x, y_rungekutt, N, h, r);

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

{

y_adams[i]=y_rungekutt[i];

}

rewit_adams (x, y_adams, N, h, r);

vivesti_rezultati (x, y_eyler, y_rungekutt, y_adams, N, f);

fclose(f);

system («pause»);

return 0;

}

Приложение 2. Программа решения краевой задачи

#include <stdio.h>

#include <iostream>

using namespace std;

const int NMAX=101;

double x[NMAX], AA[NMAX] [NMAX], BB[NMAX];

const double PI=3.14159265, EPS=1e-6;

//y''-xy'+2x=1

// у(а)=А, у(b)=B - краевые условия

int n=3, N=10;

double A=2, B=3, a=0, b=1, h=(b-a)/N;

// функции-коэффициенты диф. уравнения

// у''+p(x)*y'=f(x)

double p (double x)

{return - x;}

double f (double x)

{return 1-2*x;}

// линейная форма у''+p(x)*y'

double L (double (*f1) (double x, int k), double (*f2) (double x, int k), int k, double x)

{return f2 (x, k)+f1 (x, k)*p(x);}

// базисные функции и производные

double u (double x, int k)

{return k==0? A+(B-A)*(x-a)/(b-a): (b-x)*pow (x-a, k);}

double u1 (double x, int k)

{return k==0? (B-A)/(b-a): pow (x-a, k-1)*((b-x)*k - (x-a));}

double u2 (double x, int k)

{return k==0? 0: k*pow (x-a, k-2)*(-2*(x-a)+(k-1)*(b-x));}

// обмен местами элементов

void swap (double &a, double &b)

{double temp=a; a=b; b=temp;}

// метод гаусса, решение СЛАУ

void gauss (int n, double ans[])

{

int w[NMAX];

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

w[i]=-1;

for (int col=1, row=1; col<=n && row<=n;++col)

{

int sel=row;

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

if (fabs(AA[i] [col])>fabs (AA[sel] [col]))

sel=i;

if (fabs (AA[sel] [col])<EPS)

continue;

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

swap (AA[sel] [i], AA[row] [i]);

swap (BB[sel], BB[row]);

w[col]=row;

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

if (i!=row)

{

double c=AA[i] [col]/AA[row] [col];

for (int j=col; j<=n;++j)

AA[i] [j]-=AA[row] [j]*c;

BB[i]-=BB[row]*c;

}

++row;

}

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

ans[i]=0;

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

if (w[i]!=-1)

ans[i]=BB [w[i]]/AA [w[i]] [i];

}

// метод коллокаций (совпадений)

void Kollokatsiy (double (*u) (double x, int k), double (*u1) (double x, int k), double (*u2) (double x, int k))

{

printf («metod kollokatsiy\n»);

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

{

BB[i]=f (x[i]) - L (u1, u2,0, x[i]);

for (int j=1; j<=n;++j)

AA[i] [j]=L (u1, u2, j, x[i]);

}

printf («\n»);

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

{

for (int j=1; j<=n;++j)

printf («%7.4f», AA[i] [j]);

printf (» | %7.4f\n», BB[i]);

}

double koef[NMAX];

gauss (n, koef);

printf («resheniye SLAY\n»);

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

printf (» a[%i] =%.4f\n», i, koef[i]);

double xx=a;

printf («\n»);

printf (» X Y\n»);

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

{

double S=u (xx, 0);

for (int j=1; j<=n;++j)

S+=koef[j]*u (xx, j);

printf (» %.4f %.4f\n», xx, S);

xx+=h;

}

}

void Raznostey()

{

printf («metod kone4nix raznostey\n»);

double ck[NMAX]={}, dk[NMAX]={},

xk[NMAX]={}, yk[NMAX]={},

mk[NMAX], nk[NMAX];

xk[0]=a, xk[N]=b;

yk[0]=A, yk[N]=B;

mk[0]=-2+h*p (xk[0]),

nk[0]=1-h*p (xk[0]);

ck[0]=1/mk[0];

dk[0]=-nk[0]*A+f (xk[0])*h*h;

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

{

xk[i]=xk [i-1]+h;

mk[i]=-2+h*p (xk[i]);

nk[i]=1-h*p (xk[i]);

ck[i]=1/(mk[i] - nk[i]*ck [i-1]);

dk[i]=f (xk[i])*h*h-nk[i]*ck [i-1]*dk [i-1];

}

for (int i=N-1; i>0; - i)

yk[i]=ck [i-1]*(dk [i-1] - yk [i+1]);

printf (» X Y\n»);

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

printf (» %.4f %.4f\n», xk[i], yk[i]);

}

void qsort (double *a, int fs, int ls)

{

int l=fs, r=ls;

double f=a[(l+r)/2];

for (; l<=r;)

{

for (; a[l]<f;++l);

for (; a[r]>f; - r);

if (l<=r)

{

double t=a[l]; a[l]=a[r], a[r]=t;

++l, - r;

}

}

if (l<ls)

qsort (a, l, ls);

if (r>fs)

qsort (a, l, ls);

}

int main()

{

freopen («output.txt», «w», stdout);

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

x[i]=(b+a)/2+(b-a)/2*cos (PI*(2*i-1)/(2*n));

qsort (x, 1, n);

printf («to4ki\n»);

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

printf (» x[%i] =%.4f\n», i, x[i]);

Kollokatsiy (u, u1, u2);

Raznostey();

return 0;

}

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