1 вариант, 4 вариант, 6 вариант. Практические работы по предмету вычислительные методы

3 варианта (1 вариант, 4 вариант, 6 вариант).

ПОЛНОЕ ЗАДАНИЕ В ДЕМО ФАЙЛЕ

ЧАСТЬ ДЛЯ ПОИСКА ДУБЛИРУЮ НИЖЕ

ПРАКТИЧЕСКАЯ РАБОТА 1. ИЗМЕРЕНИЕ ФИЗИЧЕСКИХ ВЕЛИЧИН

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

1.1          Теоретическая часть

Измерения. Виды измерений

Роль измерений в науке и технике чрезвычайно велика, поскольку измерения являются основным способом получения количественной информации об окружающем нас мире и протекающих в нем процессах. Д.И.Менделеев писал: «Наука начинается с тех пор, как начинают измерять».

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

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

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

, где

,

, …

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

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

Погрешности измерений

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

Абсолютной погрешностью

-ого измерения физической величины

 называется разность

, где

- результат

-ого измерения,

- наиболее достоверное (наилучшее) значение. Согласно ГОСТ 8.207-76, за

 принимается среднее арифметическое значение всех проведенных измерений. Относительной погрешностью (точностью) 

-ого измерения называется отношение

.

Относительная погрешность часто выражается в процентах:

%.

По характеру проявления при выполнении серии измерений абсолютные погрешности делятся на  грубые, систематические и случайные. Грубые погрешности (промахи) обычно связаны с неисправностью измерительной аппаратуры, с ошибкой оператора при считывании показаний приборов, либо с резким изменением условий эксперимента, например – скачком напряжения в сети. Результаты измерений, соответствующие промахам, необходимо отбрасывать. Систематическими называются погрешности, которые при многократных измерениях остаются неизменными либо изменяются по определенному закону; они включают методические и инструментальные погрешности.

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

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

,

%.

Здесь 

- наиболее достоверное (наилучшее) значение,

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

 - относительная погрешность серии измерений. Следует отметить, что оценка величин

,

 и 

 регламентирована государственным стандартом (ГОСТ 8.207-76 «Прямые измерения с многократными наблюдениями. Методы обработки результатов наблюдений»). Для их вычислений используются соответствующие формулы теории погрешностей, которая в свою очередь базируется на теории вероятностей и математической статистике. Поскольку эти дисциплины изучаются студентами лишь в четвертом семестре, для оценки погрешностей при выполнении данной лабораторной работы используется несколько упрощенная методика.

Оценка погрешностей измерений

При обработке результатов прямых измерений выполняются следующие действия.

1.  Результаты

 измерений записываются в таблицу.

2.  Вычисляется среднее арифметическое значение результатов

 проведенных измерений

, которое принимается за наилучшее значение.

3.  Вычисляются абсолютные погрешности каждого измерения:

.

4.  Вычисляется стандартная (среднеквадратичная) погрешность:

.

Найденное значение

 сравнивается с инструментальной погрешностью прибора (

). Если

, то в качестве средней абсолютной погрешности 

 измерений принимается инструментальная погрешность:

. Если же

, то

.

5.  Вычисляется относительная погрешность серии измерений:

.

6.  Окончательный результат измерения записывается в виде

, 

.

    Обработка результатов косвенных измерений физической величины

 проводится следующим образом.

1.  Для каждой из величин

, входящих в уравнение измерения

, проводится обработка в соответствии с перечисленными выше п.1…4 (иначе говоря, для каждой из величин  

 находятся наилучшие значения и средние абсолютные погрешности).

2.  Вычисляется наилучшее значение измеряемой величины

:

(здесь

- наилучшие значения величин

).

3.  Вычисляется относительная погрешность серии измерений по формуле

.

Здесь символическая запись

 обозначает частную производную по переменной

 прологарифмированного уравнения измерения (примеры вычисления частных производных рассматриваются в отчете по лабораторной работе «Измерение плотности прямоугольного параллелепипеда», включенного в «Приложение» в качестве образца). Численное значение частной производной по переменной

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

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

 переменных.

4.   Вычисляется средняя абсолютная погрешность величины

:

.

5.  Окончательный результат измерения записывается в виде

, 

.

1.3 Порядок выполнения

Задача: Плотность твердого тела находится путем косвенного измерения, т.е. численные значения его массы и объема подставляются в формулу

.

Здесь

 - масса тела (предполагается известной),

 - объем тела как функция его

 линейных размеров.

Решение

1.             Изучите правила приближенных вычислений (Приложение А).

2.             Изучите отчет к практической  работе «Измерение плотности прямоугольного параллелепипеда», приведенный в «Приложения Б» в качестве образца.

3.             Изучите методику измерения линейных размеров тела штангенциркулем.

4.             Сделайте эскиз тела, предложенного преподавателем, с обозначением всех величин, подлежащих измерению.

5.              Составьте таблицу результатов измерений.

6.             Проведите необходимые измерения не менее трех раз и результаты поместите в таблицу.

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

8.              

 кг/м3, 

%.

ПРИЛОЖЕНИЕ А

П.1. Правила вычислений при обработке результатов измерений

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

           Значащими цифрами числа называются все отличные от нуля цифры и нули, расположенные между этими цифрами и в конце числа. Например, числа 0,072 и 8,4 имеют по две значащие цифры, число 6,40 – три значащие цифры. При считывании показаний измерительного прибора цифра числа является верной,   если половина цены разряда, в котором она находится, не меньше инструментальной погрешности. Например, при измерении длины школьной линейкой получен результат 23,1 см. Последняя цифра верна, поскольку она находится в разряде миллиметров (половина цены этого разряда составляет 0,5 мм, что не меньше инструментальной погрешности линейки).

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

1.    Если первая отбрасываемая цифра больше 4, то последняя сохраняемая цифра увеличивается на единицу. Например, округляя число 27,3763 до сотых, следует записать 27,38.

2.    Если первая отбрасываемая цифра равна 4 или меньше 4, то последняя сохраняемая цифра не изменяется. Например, округляя число 13847 до сотен, записывают 138∙102.

3.    Если отбрасываемая часть числа состоит из одной цифры 5, то число округляют так, чтобы последняя сохраняемая цифра была четной. Например, при округлении до десятых 23,65≈23,6, но 17,75≈17,8.

Производя различные математические действия над приближенными числами, руководствуются следующими правилами.

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

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

3.    При возведении в степень (целую или дробную) в результате оставляют столько значащих цифр, сколько их в основании степени.

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

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

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

ПРИЛОЖЕНИЕ Б

«Измерение плотности прямоугольного параллелепипеда»

           Приборы и принадлежности: штангенциркуль, весы лабораторные.

           Эскиз тела:

Уравнение измерения:

,                                                                  (1.1)

где 

- плотность, 

 - масса,

 - линейные размеры тела.

           Результаты измерений:

№ изм.

, мм

, мм

, мм

,  г

1

38,9

29,4

17,5

160,3

2

39,0

30,1

17,4

160,3

3

39,1

29,9

17,5

160,3

Вычисление наилучших значений линейных размеров:

мм,

мм,

мм.

           Вычисление абсолютных погрешностей:

мм, аналогично находим

,

мм;

мм, аналогично находим

мм,

мм;

, аналогично находим

мм, 

.

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

мм; поскольку

больше инструментальной погрешности штангенциркуля (0,05 мм),

мм;

мм; поскольку

, 

мм;

мм; поскольку

, 

мм.

Учитывая, что масса тела предполагается известной, считаем, что

Вычисление наилучшего значения плотности тела:

кг/м

3

.

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

 уравнения измерения (1.1), прологарифмированного по натуральному основанию:

.

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

 по переменной

, остальные переменные полагаем константами и находим:

.

Аналогично имеем:

, 

, 

.

В соответствии с этим получим формулу для вычисления относительной погрешности:

.

В результате подстановки численных значений имеем (с учетом правил приближенных вычислений):

что составляет 0,7%.

Вычисляем среднюю абсолютную погрешность:

 кг/м

3

.         

           Таким образом,

кг/м

3

,

.

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

Практическое задание № 4. Решение задач систем линейных уравнений

1 Цель занятия

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

- Изучить итерационные методы решения СЛАУ (метод простых итераций, метод Зейделя).

- Составить алгоритм решения СЛАУ указанными методами, применимый для организации вычислений на ЭВМ.

- Составить программу решения СЛАУ по разработанному алгоритму.

- Численно решить тестовые примеры и проверить правильность работы программы. Сравнить трудоемкость решения методом простых итераций и методом Зейделя.

2 Задание

Вычислить с точностью ε=10-3 СЛУодним из методов (таблица 2.1).

Таблица .1

Вариант

СЛУ

Метод

1

4,5x1-1,7x2+3,5x3=23,1

x1+2,3x2-1,1x3=11,8

x1+2,5x2+4,7x3=4

Зейделя

2

0,23x1-0,04x2+0,21x3=1,24

0,5x1+0,02x2-3,1x3=12,4

2x1+0,1x2-5x3=10,7

простой итерации

3

5x1-7x2+5x3=3,1

2x1+3x2-0,1x3=1,8

3x1+5x2+7x3=14

Зейделя

4

4x1-x2+5x3=20,1

x1+2,3x2-1,1x3=0,8

0,1x1+0,5x2+0,7x3=0,4

простой итерации

5

0,5x1-0,7x2+0,5x3=0,1

0,1x1+0,3x2-0,1x3=0,8

0,1x1+0,5x2+0,7x3=0,7

Зейделя

6

x1-7x2+5x3=20

2x1+3x2-x3=8

x1+8,5x2+4x3=24

простой итерации

7

0,6x1-0,5x2+0,4x3=0,3

0,5x1+0,2x2-0,9x3=5

0,8x1+1,1x2+2,3x3=0,5

Зейделя

8

0,1x1-1,3x2+1,5x3=7,1

1,2x1+1,4x2-1,1x3=0,7

2,1x1+2,5x2+2,7x3=7,7

простой итерации

9

0,3x1-0,3x2+5,2x3=2

0,2x1+0,3x2-0,1x3=3

5x1+6x2+7x3=4

Зейделя

10

6,2x1-7,1x2+4x3=30,7

0,2x1+0,4x2-1,1x3=11

8x1+5x2+3x3=14

простой итерации

3 Указания к выполнению заданий

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

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

Для решения СЛАУ методом простых итераций преобразуем систему от первоначальной формы Ax = b или

a11x1 + a12 x2 + … + a1nxn = b1 ,

a21x + a22x2 + … + a2nxn = b2 ,                                                           (1)

an1x1 + an2x2 + … + annxn=bn

квиду  x = Bx + c.                                                                               (2)

где B – квадратная матрица с элементами bij (i, j = 1, 2, …, n),

     c – вектор-столбец с элементами ci (i = 1, 2, …, n).

В развернутой форме записи система (2) имеет следующий вид:

x1 = b11x1 + b12x2 + b13x3 + … + b1nxn + c1

x2 = b21x1 + b22x2 + b23x3 + … + b2nxn + c2

. . . . . . . . . . . . . . . . .

xn = bn1x1 + bn2x2 + bn3x3 + … + bnnxn + cn

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

Можно, например, преобразовать систему (1) следующим образом

x1 = (b1 - a11x1 – a12x2 - … - a1nxn)/a11 + x1,

x2 = (b2 - a21x1 – a22x2 - … - a2nxn)/a22 + x2,

 …..

xn = (bn - an1x1 – an2x2 - … - annxn)/ann + xn

если диагональные элементы матрицы А отличны от нуля.

Можно преобразовать систему (1) в эквивалентную ей систему

x = (E-A)x+b

Задав произвольным образом столбец начальных приближений х0 = (x10 , x20 , … , xn0 ), подставим их в правые части системы (2) и вычислим новые приближения x1 = (x11 , x21 , … , xn1 ), которые опять подставим в систему (2) и т.д. Таким образом, организуется итерационный процесс

xk=Bxk-1+c, k=1, 2, …. ,.

Известно, что система (1) имеет единственное решение х* и последовательность {xk} сходится к этому решению со скоростью геометрической прогрессии, если ||B||<1 в любой матричной норме.

Метод Зейделя

Метод Зейделя является модификацией метода простых итераций. Суть его состоит в том, что при вычислении следующего элемента в формуле xk=Bxk-1+c, k=1,2, …. используются вместо xi(k-1), уже вычисленные ранее значение xi(k).

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

Схема алгоритма аналогична схеме метода простых итераций.

Самый простой способ приведения системы к виду, удобному для итераций, состоит в следующем. Из первого уравнения системы выразим неизвестное x1:

x1 = a11–1(b1 – a12x2 – a13x3 – … – a1nxn),

из второго уравнения – неизвестное x2:

x2 = a21–1(b2 – a22x2 – a23x3 – … – a2nxn),

и т. д. В результате получим систему

x1 = b12x2 + b13x3 + … + b1,n–1xn–1 + b1nxn+ c1,

x2 = b21x1 + b23x3 + … + b2,n–1xn–1 + b2nxn+ c2 ,

x3 = b31x1 + b32x2 + … + b3,n–1xn–1 + b3nxn+ c3 ,

. . . . . . . . . . . . . . . . . . . . . . .

xn = bn1x1 + bn2x2 + bn3x3 + … + bn,n–1xn–1 + cn ,

в которой на главной диагонали матрицы B находятся нулевые элементы. Остальные элементы выражаются по формулам

bij = –aij/aii, ci = bi/aii (i, j = 1, 2, …, n, j ≠ i)    (2.4)

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

Введем нижнюю B1 (получается из B заменой нулями элементов стоявших на главной диагонали и выше ее) и верхнюю B2 (получается из B заменой нулями элементов стоявших на главной диагонали и ниже ее) треугольные матрицы.

Заметим, что B=B1+B2 и поэтому решение x исходной системы удовлетворяет равенству

                                   x = B1x + B2x + c                                          (5)

Выберем начальное приближение x(0)=[x1(0), x2(0), …, xn(0)]. Подставляя его в правую часть равенства при верхней треугольной матрице B2 и вычисляя полученное выражение, находим первое приближение

x(1) = B1x(0) + B2x(1)  (6)

Подставляя приближение x(1), получим

x(2) = B1x(1) + B2x(2)  (7)

Продолжая этот процесс далее, получим последовательность x(0), x(1), …, x(n), … приближений к вычисляемых по формуле

                                     x(k+1) = B1(k+1) + B2(k)+c                          (8)

или в развернутой форме записи

x1(k+1) = b12x2(k) + b13x2(k) + … + b1nxn(k) + c1,

x2(k+1) = b21x1(k+1) + b23x3(k) + … + b2nxn(k) + c2,

x3(k+1) = b31x1(k+1) + b32x2(k+1) + … + b3nxn(k) + c3 ,

. . . . . . . . . . . . . . . . . . . . . . . . . .

xn(k+1) = bn1x1(k+1) + bn2x2(k+1) + bn3x3(k+1) + … + cn .

Объединив приведение системы к виду, удобному для итераций и метод Зейделя в одну формулу, получим

xi(k+1) = xi(k) – aii–1(∑j=1i–1 aijxj(k+1) + ∑j=1n aijxi(k) – bi) (2.9)

Тогда достаточным условием сходимоти метода Зейделя будет условие доминированния диагональных элементов в строках или столбцах матрицы A, т.е.

aii>ai1+ ……..+ain для всех i=1,…n,

или ajj>a1j+…….+anj для всех j=1,….,n.

Методы простой итерации и Зейделя сходятся примерно так же, как геометрическая прогрессия со знаменателем ||B||

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

function [converge, X, iterations, closure] = SimpleIterationsMethod( A, b, eps, max )

%Solves the system of linear equations

B = zeros(size(A,1), size(A, 2));

c = b;

for i = 1 : size(A, 1)

c(i) = c(i) / A(i,i);

for j = 1 : size(A, 2)

if (i == j)

B(i, j) = 0;

else

B(i, j) = - A(i, j) / A(i, i);

end;

end;

end;

converge = true;

k = 1;

x0 = c;

xk = B * x0 + c;

norm_no = 1;

coef = 0;

if (norm(B, norm_no) > 1)

norm_no = 2;

if (norm(B, norm_no) > 1)

norm_no = Inf;

if (norm(B, norm_no) > 1)

norm_no = 1;

coef = 1;

end;

end;

end;

if (coef ~= 1)

coef = norm(B, norm_no) / (1 - norm(B, norm_no)) ;

end;

while (norm(xk - x0, norm_no) * coef >= eps)

x0 = xk;

xk = B * x0 + c;

k = k + 1;

if (k == max && coef == 1)

converge = false;

break;

end;

end

X = xk;

iterations = k;

closure = A * X - b;

end

function [ x, best_k, best_t, closure, converge ] = SimpleIterationsMethod2( A, b, eps, max_iterations, t_min, t_max)

%Solves the system of linear equations

N = size(A, 1);

best_t = 0;

best_k =Inf;

step = (t_max - t_min) / 1000;

fileId = fopen('report.txt', 'w');

for t = t_min : step : t_max

if (t == 0)

continue;

end;

B = eye(N) - t * A;

c = t * b;

norm_no = 1;

coef = 0;

if (norm(B, norm_no) > 1)

norm_no = 2;

if (norm(B, norm_no) > 1)

norm_no = Inf;

if (norm(B, norm_no) > 1)

norm_no = 1;

coef = 1;

end;

end;

end;

if (coef ~= 1)

coef = norm(B, norm_no) / (1 - norm(B, norm_no)) ;

end;

k = 1;

x0 = c;

xk = B * x0 + c;

converge = true;

while (norm(xk - x0, norm_no) * coef >= eps )

x0 = xk;

xk = B * x0 + c;

k = k + 1;

for i = 1 : size(xk)

if (isnan(xk(i)) || isinf(xk(i)))

converge = false;

break;

end;

end;

if (converge == false || (coef == 1 && k == max_iterations))

converge = false;

break;

end;

end;

fprintf(fileId, '\tX when t = %.4f, k = %d\n', t, k);

if (converge)

fprintf(fileId, '%.4f\n', xk);

else

fprintf(fileId, 'Converge == false');

end;

fprintf(fileId, '\n\n');

if (k < best_k && converge == true)

best_k = k;

best_t = t;

x = xk;

end;

end;

converge = best_k < Inf;

if (converge)

closure = A * x - b;

else

closure = Inf;

x = Inf;

end;

fclose(fileId);

end

function [converge, X, iterations, U] = ZeidelMethod( A, b, eps, max )

B = zeros(size(A,1), size(A, 2));

c = b;

for i = 1 : size(A, 1)

c(i) = c(i) / A(i,i);

for j = 1 : size(A, 2)

if (i == j)

B(i, j) = 0;

else

B(i, j) = - A(i, j) / A(i, i);

end;

end;

end;

useMax = false;

norm_no = 1;

if (norm(B, norm_no) > 1)

norm_no = Inf;

if (norm(B, norm_no) > 1)

useMax = true;

end;

end;

n = norm(B, norm_no);

x = c;

k = 0;

while (true)

for i = 1 : size(A, 1)

var = 0;

for j = 1 : size(A, 2)

if (j ~= i)

var = var + B(i, j) * x(j);

end;

end;

x(i) = (c(i) + var) ;

end;

k = k + 1;

if (norm(A*x - b, norm_no) <= eps)

break;

end;

if (useMax && k >= max)

break;

end;

end;

X = x;

iterations = k;

if (useMax && k == max)

converge = false;

U = Inf;

else

converge = true;

U = A*X - b;

end;

end

4 Контрольные вопросы

1. Опишите суть методы простой итерации для решения СЛУ

2. В чем отличие и сходство при решении СЛУ методом Зейделя

3. Какие начальные условия необходимо принять для решения СУ методом простой итерации

4. Опишите алгоритм решения методом Зейделя

5 Выводы по заданию

Практическая работа 6. Основы решения дифференциальных уравнений

Цель занятия:

- закрепить усвоение теоретического материала по данной теме через решение упражнений;

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

6.1 Задание

Решить задачу Коши для дифференциального уравнения y’=f(x,y) на отрезке [a;b] при заданном начальном условии y(a)=y0 и шаге интегрирования h методом Эйлера с применением «ручных» вычислений с шагом 2h. Варианты задания представлены в таблице 6.1. Пользуясь таблицей 6.1, сделать ручную прикидку графика интегральной кривой на бумаге.

Таблица 6.1 –Варианты заданий

6.2 Теоретические сведения

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

y’=f(x,y) (6.1)

Эта задача известна, как задача Коши: найти решение уравнения (6.1) в виде функции y(x), удовлетворяющей начальному условию

                               y(x0) = y0                                                    (6.2)

Геометрически это означает, что требуется найти интегральную кривую y=y(x), проходящую через заданную точку М0 (x0,y0), при выполнении равенства (6.1).

Существует несколько классов дифференциальных уравнений 1-го порядка, для которых решение может быть найдено аналитически. Но даже для таких уравнений решение не всегда удается довести до вида y=y(x). Многие же дифференциальные уравнения, к которым приводят математические модели реальных процессов, не могут быть решены аналитически.

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

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

2) графические методы, дающие приближенное решение в виде графика;

3) численные методы, когда искомая функция получается в виде таблицы.

Метод Эйлера

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

Пусть дано уравнение (6.1) с начальным условием (6.2), т.е. поставлена раздача Коши. Вначале найдем простейшим способом приближенное значение решения в некоторой точке x1= x0+ h ,где h – достаточно малый шаг.

Заметим, что уравнение (6.1) совместно с начальным условием (6.2) задают направление касательной к искомой интегральной кривой в точке М0(x0,y0). Двигаясь вдоль этой касательной, получим приближенное значение решения в точке

                             y1= y0+hf x0;y0)                                          (6.3)

Аналогично, найдем приближенное значение решения в точке x2= x1+h, и т.д.

Продолжая эту идею, построим систему равностоящих точек xi=x0+ih, i=0,..,n.

Получение таблицы значений искомой функции y(x) по методу Эйлера заключается в циклическом применении пары формул:

(6.4)

Геометрическая иллюстрация метода Эйлера:

Рисунок 6.1 – Метод Эйлера

Вместо кривой в реальности получается совокупность прямых – ломаная Эйлера. Методы численного интегрирования дифференциальных уравнений, в которых решение получается от одного узла к другому, называются пошаговыми. Метод Эйлера – простейший пошаговый метод. Отметим, что оценка погрешности метода при таком элементарном рассмотрении невозможна даже на первом шаге. Кроме того, особенностью любого пошагового метода является то, что, начиная со второго шага, исходное значение yi в формуле (6.4) само является приближенным, т.е. погрешность на каждом шаге систематически возрастает. Наиболее используемым методом оценки точности, как метода Эйлера, так и других пошаговых методов приближенного численного интегрирования обыкновенных дифференциальных уравнений является способ двойного прохождения заданного отрезка с шагом h и с шагом h/2. Совпадение соответствующих десятичных знаков в полученных двумя способами результатах дает основание считать их верными.

Пример:

Решить методом Эйлера дифференциальное уравнение y¢=cos(y)+ 3x c начальным условием y(0) = 1,3 на отрезке [0;1] применив h=0,2. Решение: Имеем f(x, y)=cos(y)+ 3x . Составим таблицу значений функции f(x,y) с шагом h и h/2.

Таблица 6.1

При составлении таблицы проводились следующие вычисления:

Если h=0,2:

1) х0=0, у0=1,3 из начального условия;

2) х1=0,1,

3) x2=0,2

и т.д.

Аналогичные вычисления проводились и для h=0,1.

Таким образом, приближенное решение уравнения получаем в виде таблицы. Построим ломаную Эйлера для h=0,2 и h=0,1 в одной системе координат.

Рисунок 6.2 – Метод Эйлера на примере

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

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

Альтернативный путь открывает метод Рунге-Кутта, названный по имени его создателей.

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

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

                        (6.5)

                                 (6.6)

Общий недостаток методов Рунге-Кутта – отсутствие простых способов оценки погрешности метода. Погрешность на одном шаге оценить сравнительно не трудно, гораздо труднее оценить накопление погрешностей на протяжении многих шагов. Широко используемый на практике для этих методов способ контроля точности – двойной счет: вычисляем решение дифференциального уравнение с шагом h и h/2 , а потом сравниваем полученные результаты.

Пример:

Решить дифференциальное уравнение y¢= y(1- x) на отрезке [0;0.5] с начальным условием у(0)=1 и шагом h=0.05.

Решение:

Сначала решим это уравнение аналитически: y’= y (1- x) - уравнение с разделяющимися переменными.

Применим начальное условие y(0) =1, получим:

Таким образом, частное решение данного уравнения, удовлетворяющее заданному начальному условию:

Пользуясь этой формулой, можно получить таблицу «точного» решение уравнения. Найдем приближенное решение дифференциальное уравнение по методу Рунге-Кутта. Проведем последовательные вычисления по формулам (6.5), (6.6):

Имеем:

Тогда

Подставим найденные значения в формулу (6.5):

Поскольку вычисления достаточно громоздки и трудоемки, то численные решения заданного уравнения можно найти с помощью программы на компьютере. Для сравнения результатов построим таблицу, в которой укажем численные решения, полученные по методу Эйлера, методу Рунге-Кутта и «точное решение».

Из таблицы видно, что результаты, получения по методу Рунге-Кутта практически совпадают с «точным» решением уравнения, в отличие от соответствующих значений, полученных по методу Эйлера.

6.3 Контрольные вопросы

1. Что является решением дифференциального уравнения?

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

3. В какой форме получается приближенное решение дифференциального уравнения по методу Эйлера?

4. В чем основная идея метода Рунге-Кутта?

5. В чем отличие одношаговых методов Эйлера и Рунге-Кутта?