Численные методы. Номер 1.. 1. Приближенное вычисление определенных интегралов приближенное вычисление интеграла по формулам трапеций и Симпсона, оценка погрешностей вычисления Формула трапеций имеет вид
 Скачать 4.47 Mb.
|
|
y 0 z 0 1 0 ) ), ( , ( x z x y x f 56 y xz x z x y x g 09 0 ) ), ( , ( y 0 =0.7 z 0 =1/11 2 0 ) ), ( , ( x z x y x f 3 x e x z x y x g 4 0 07 0 ) ), ( , ( y 0 =1.3 z 0 =0.125 3 0 ) ), ( , ( x z x y x f 8x y x xy x z x y x g 2 11 0 ) ), ( , ( y 0 =-1 z 0 =1/11 4 ) ), ( , ( x x z x y x f y xz x z x y x g 1 0 ) ), ( , ( y 0 =0.7 z 0 =1/12 5 y x z x y x f ) ), ( , ( x e x z x y x g 4 0 08 0 ) ), ( , ( y 0 =0.3 z 0 =1/9 6 ) ), ( , ( y x x z x y x f y x xz x z x y x g 2 12 0 ) ), ( , ( y 0 =1.6 z 0 =1/12 7 5 0 ) ), ( , ( x z x y x f x e x yz x y x g 4 0 09 0 ) ), ( , ( y 0 =0 z 0 =0.1 8 y x z x y x f ) ), ( , ( x e x z x y x g 5 0 1 0 ) ), ( , ( y 0 =1.3 z 0 =0.1 9 х y x xz x z x y x g 2 13 0 ) ), ( , ( y 0 =0.8 z 0 =1/13 10 z y x z x y x f 2 ) ), ( , ( z y x z x y x g 2 ) ), ( , ( y 0 =-1 Задача опадении тел Задание 13. Тело массой 5 кг начинает падать с высоты Н, значит 50 м, 0 0 t v . Вычислить время падения тела до земли с учетом сопротивления воздуха, пропорционального скорости падения, и без учета сопротивления воздуха. Решение. Согласно второму закону Ньютона Закон сопротивления воздуха пропорционален скорости v k F d , k=0.0001, начальные условия 0 ) ( 0 t t v и g m v k mg t a t t 0 Ищем численное решение Зная начальные условия приищем решение в точке t t t 0 1 (6.5) t g t t a t v t v 0 ) ( ) ( ) ( 0 0 1 (6.6) t t v t y t y ) ( ) ( ) ( 1 0 1 (6.7) m t v k mg t a ) ( ) ( 1 1 (6.8) Затем вычисляем решение в точке t t t 1 Решение в Microsoft Excel: первый столбец таблицы (А) отводим для значений времени, в которые вычисляем скорость и высоту тела, во втором столбце будем вычислять скорость тела, в третьем – высоту, в четвертом ускорение. В пятом и шестом столбцах записываем исходные данные задачи. В третьей строке записываем начальные условия задачи отсчет времени начинается с нуля, начальная скорость равняется нулю, начальная высота, с которой падает телом, начальное ускорение равняется 9,81. В строках первого столбца записываем значения времени, в которых ищем решение задачи. Выбираем шаг изменения времени 0,1 до значения 3,3. В столбце В записываем формулу (6.6) для вычисления скорости, в столбце С формулу (6.7) для вычисления высоты, в четвертом – формулу (6.8) для вычисления ускорения (рис. 6.6, 6.8). Одновременно выполняем расчет времени падения тела без учета сопротивления воздуха. Рис. 6.6. Фрагмент таблицы с решением задачи опадении тела в режиме отображения формул 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 с сопротивлением воздуха t v y a дельта t 0,1 0 0 50 нач высота 50 масса к 0,3 =B5+D5*0,1 =C5-B6*0,1 =($F$4*9,81-$F$5*B6)/$F$4 0,4 =B6+D6*0,1 =C6-B7*0,1 =($F$4*9,81-$F$5*B7)/$F$4 0,5 =B7+D7*0,1 =C7-B8*0,1 =($F$4*9,81-$F$5*B8)/$F$4 0,6 =B8+D8*0,1 =C8-B9*0,1 =($F$4*9,81-$F$5*B9)/$F$4 0,7 =B9+D9*0,1 =C9-B10*0,1 =($F$4*9,81-$F$5*B10)/$F$4 0,8 =B10+D10*0,1 =C10-B11*0,1 =($F$4*9,81-$F$5*B11)/$F$4 0,9 =B11+D11*0,1 =C11-B12*0,1 =($F$4*9,81-$F$5*B12)/$F$4 1 =B12+D12*0,1 =C12-B13*0,1 =($F$4*9,81-$F$5*B13)/$F$4 1,1 =B13+D13*0,1 =C13-B14*0,1 =($F$4*9,81-$F$5*B14)/$F$4 1,2 =B14+D14*0,1 =C14-B15*0,1 =($F$4*9,81-$F$5*B15)/$F$4 1,3 =B15+D15*0,1 =C15-B16*0,1 =($F$4*9,81-$F$5*B16)/$F$4 1,4 =B16+D16*0,1 =C16-B17*0,1 =($F$4*9,81-$F$5*B17)/$F$4 1,5 =B17+D17*0,1 =C17-B18*0,1 =($F$4*9,81-$F$5*B18)/$F$4 1,6 =B18+D18*0,1 =C18-B19*0,1 =($F$4*9,81-$F$5*B19)/$F$4 1,7 =B19+D19*0,1 =C19-B20*0,1 =($F$4*9,81-$F$5*B20)/$F$4 1,8 =B20+D20*0,1 =C20-B21*0,1 =($F$4*9,81-$F$5*B21)/$F$4 1,9 =B21+D21*0,1 =C21-B22*0,1 =($F$4*9,81-$F$5*B22)/$F$4 2 =B22+D22*0,1 =C22-B23*0,1 =($F$4*9,81-$F$5*B23)/$F$4 2,1 =B23+D23*0,1 =C23-B24*0,1 =($F$4*9,81-$F$5*B24)/$F$4 2,2 =B24+D24*0,1 =C24-B25*0,1 =($F$4*9,81-$F$5*B25)/$F$4 2,3 =B25+D25*0,1 =C25-B26*0,1 =($F$4*9,81-$F$5*B26)/$F$4 2,4 =B26+D26*0,1 =C26-B27*0,1 =($F$4*9,81-$F$5*B27)/$F$4 2,5 =B27+D27*0,1 =C27-B28*0,1 =($F$4*9,81-$F$5*B28)/$F$4 2,6 =B28+D28*0,1 =C28-B29*0,1 =($F$4*9,81-$F$5*B29)/$F$4 2,7 =B29+D29*0,1 =C29-B30*0,1 =($F$4*9,81-$F$5*B30)/$F$4 2,8 =B30+D30*0,1 =C30-B31*0,1 =($F$4*9,81-$F$5*B31)/$F$4 2,9 =B31+D31*0,1 =C31-B32*0,1 =($F$4*9,81-$F$5*B32)/$F$4 3 =B32+D32*0,1 =C32-B33*0,1 =($F$4*9,81-$F$5*B33)/$F$4 3,1 =B33+D33*0,1 =C33-B34*0,1 =($F$4*9,81-$F$5*B34)/$F$4 3,2 =B34+D34*0,1 =C34-B35*0,1 =($F$4*9,81-$F$5*B35)/$F$4 3,3 =B35+D35*0,1 =C35-B36*0,1 =($F$4*9,81-$F$5*B36)/$F$4 Рис. 6.7. Решение задачи опадении тела с учетом сопротивления воздуха Проанализируем полученный результат. Положение тела над землей (столбец С ) убывает, так как тело падает с некоторой высоты. Через примерно 3,1 сек. тело упадет на землю. С этого момента формулы определения высоты дают уже отрицательное значение. Это значит, что вычисления уже не имеют смысла. Сравнивая результаты, можем утверждать, что учет сопротивления воздуха при падении тела вносит незначительное изменение в результат. Это говорит о том, что математическая модель падающего тела без учета сопротивления воздуха хорошо описывает его падение, и учет сопротивления воздуха не вносит значительных изменений. Решение в пакете MathCAD Для решения обыкновенных дифференциальных уравнений применим функцию odesolve, которая применяется для решения уравнений любого порядка. В связи стем, что ось ординат обычно направлена вверх, а падение идет вниз в уравнении перед второй производной (ускорением) ставится минус. Решение приведено на рис. 6.8. Рис. 6.8. Решение в MathCAD задачи опадении тела с учетом сопротивления воздуха ВАРИАНТЫ ЗАДАНИЯ ПО ТЕМЕ 6 Тело массой m начинает падать с высоты Н. Вычислить время падения тела до земли с учетом сопротивления воздуха пропорциональным степени скорости падения с d v k F Решение выполнить в табличном процессоре Microsoft Excel и пакете математических расчетов MathCAD. В отчете привести систему дифференциальных уравнений, описывающих падение тела формулы Эйлера для вычисления численного решения системы дифференциальных уравнений первого порядка численное решение системы уравнений по формулам Эйлера в Microsoft Excel (в режимах отображения чисел и формул графическое представление решения приближенное решение, полученное в пакете Mathcad с использованием функции odesolve. Вариант 1. m=10 кг, Нм. Вариант 2. m=50 кг, Нм. Вариант 3. m=75 кг, Нм. Вариант 4. m=25 кг, Нм. Вариант 5. m=100 кг, Нм. Вариант 6. m=10 кг, Нм. Вариант 7. m=10 кг, Нм. Вариант 8. m=10 кг, Нм. Вариант 9. m=10 кг, Нм. Вариант 10. m=10 кг, Нм. ТЕМА 7. ЧИСЛЕННОЕ РЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ВТОРОГО ПОРЯДКА. ЗАДАЧА КОШИ Рассмотрим обыкновенные дифференциальные уравнения второго порядка x y x y x f x y , , на отрезке изменения аргумента x [a,b] с начальными условиями y(a)=y a , y'(a)=y' a . Они содержат вторую производную неизвестной функции. В решении дифференциального уравнения второго порядка присутствуют две неизвестные постоянные, которые определяются изначальных условий С помощью подстановки обыкновенное дифференциальное уравнение второго порядка сводится к системе двух обыкновенных дифференциальных уравнений первого порядка 1 1 1 , , y y x G y y y . К каждому из этих уравнений применим метод Эйлера или метод Рунге-Кутта. Задание 14. Задано обыкновенное дифференциальное уравнение второго порядка с начальными условиями и 5 , 0 ) 1 ( y . Найти приближенное решение уравнения на промежутке x [1, 2]. С помощью подстановки z y z y , заданное обыкновенное дифференциальное уравнение второго порядка сводится к системе двух обыкновенных дифференциальных уравнений первого порядка x z y x z z y sin с начальными условиями 5 , 0 ) 1 ( y и 5 , 0 ) 1 ( z . Применяя к каждому из уравнений системы метод Эйлера, можем записать систему в следующем виде i i i i i i i i i x z y x h z z z h y y sin 1 1 , (7.1), где x 0 =1; y 0 =0,5; z 0 =-0,5. Вычисления производим последовательно по каждой из формул системы, изменяя значение i от нуля до n, где n – число точек разбиения отрезка, на котором определяется решение. Каждое следующее значение аргумента вычисляется по формуле h x x i i 1 , где h – шаг изменения аргумента. Решение в табличном процессоре Microsoft Excel будем проводить в следующей последовательности надписываем столбцы A, B, C как x, y, y’. В ячейку Е записываем значение шага изменения аргумента. Заполняем ячейки столбца А значениями аргументах от начального значения 1 до конечного 2 шагом 0,1 (рис. 7.1). Набранные формулы нужно скопировать до конечного значения аргумента - значениях. В результате копирования Microsoft Excel выведет полученные значения (рис. 7.2). Рис. 7.1. Начальный этап решения дифференциального уравнения второго порядка Рис. 7.2. Ввод формул для вычисления y, y Во второй строке подзаголовками записываем заданные начальные значения. В третьей строке записываем формулы для вычисления y, y’ по формулам системы (7.1) рис. 7.3). При записи формул нужно не забывать ставить абсолютный адрес ячейки, содержащей значение шага. Полученное решение нужно проиллюстрировать графически. Для этого выделяем диапазоны ячеек, содержащие значения аргумента, функции и первой производной, вызываем мастер диаграмм и строим диаграмму типа точечная (рис. 7.4). Рис. 7.3. Вывод результата вычислений Рис. 7.4. Графическое решение уравнения Для решения задачи в пакете MathCAD применяем функцию odesolve (ordinary differentional equation solve). Эта функция может искать решения обыкновенных дифференциальных уравнений любого порядка. Функция имеет три аргумента имя аргумента искомой функции, правый конец интервала интегрирования, число шагов разбиения интервала интегрирования. Результатом работы функции является вычисленная функция. Решение приведено на рис. 7.5. Рис. 7.5. Графическое решение уравнения с использованием функции odesolve Решение дифференциального уравнения может быть получено и с использованием функции решения обыкновенного дифференциального уравнения rkfixed. Результатом работы функции является матрица решений системы дифференциальных уравнений методом Рунге-Кутта с постоянным шагом. Пакет выводит матрицу с решением, в первом столбце которой содержатся значениях, для которых вычислено решение, во втором столбце – значения искомой функции в заданных точках, в третьем – значения первых производных в тех же точках. Для получения графической формы решения вызывается палитра Graph пакета, активизируется кнопка плоского графика. В маркере у оси абсцисс задается имя нулевого столбца матрицы решения (в котором содержатся значения аргументах, в маркере у оси ординату задается имя первого столбца матрицы решения, содержащего полученные значения искомой функции. Для задания номера столбца удобно воспользоваться значком указания столбца палитры Matrix <> . На рис. 7.6 приведено описанное решение. Рис. 7.6. Графическое решение уравнения с использованием функции rkfixed Решения, полученные всеми способами в обоих пакетах, совпадают. Чтобы проверить точность полученного решения, найдем точное решение уравнения. Точное решение имеют линейные дифференциальные уравнения второго порядка с постоянными коэффициентами, те. уравнения вида 0 2 1 y a y a y , где а, а – постоянные величины. Общим решением уравнения является функция С 1 2 1 ) ( , где С, С – постоянные, k 1 , k 2 – корни характеристического уравнения 0 2 1 2 a k a k . (Если корни равные, то решение имеет вид С 1 2 1 ) ( , если корней вещественных нет, те. i k 2 , 1 , то решение имеет вид x C x C e x y x sin cos 2 1 ). Задание 15. Найти решение обыкновенного дифференциального линейного уравнения второго порядка с постоянными коэффициентами 0 3 y y y на диапазоне изменения аргумента x [-0.5, 1] при начальных условиях y(-0.5)=-1, y′(-0.5)=0.5. Решение. Составим характеристическое уравнение 0 1 3 2 k k . Корни уравнения вычисляются по формуле 1 4 9 2 3 2 , 1 k . Корни различны по величине, следовательно, решение имеет вид x x e C e C x y 2 5 2 3 2 2 5 2 3 1 ) ( . Для определения постоянных С, С применим начальные условия уравнения 5 0 2 5 2 3 2 5 0 2 5 2 3 1 1 e C e C (8.4). Вычислим первую производную уравнения и подставим ее во второе начальное условие x x e C e C x y 2 5 2 3 2 2 5 2 3 1 2 5 2 3 2 5 2 3 ) ( , следовательно, 5 0 2 5 2 3 2 5 0 2 5 2 3 1 2 5 2 3 2 5 2 3 5 0 e C e C (7.5). Те. для определения постоянных нужно решить систему линейных алгебраических уравнений (8.4) и (8.5). Решение выполним в пакете MathCAD. Решение приведено на рис. 7.7. Сопоставим полученное точное решение с приближенным, вычисленным с использование функций odesolve (рис. 7.8) ирис Рис. 7.7. Точное решение уравнения Рис. 7.8. Решение уравнения с использованием функции odesolve Рис. 7.9. Графическое решение уравнения с использованием функции rkfixed ВАРИАНТЫ ЗАДАНИЙ ПО ТЕМЕ 7 ЗАДАНИЕ 13. Задано обыкновенное дифференциальное уравнение второго порядка ) , , ( x y y f y с начальными условиями 0 ) ( y a y и 0 ) ( y a y . Найти точное решение уравнения. Найти численное решение уравнения методом Эйлера на промежутке x [a, b]. Оценить максимальную погрешность вычисления. Решение выполнить в Microsoft Excel и пакете Mathcad. В отчете привести точное решение уравнения формулу Эйлера вычисления численного решения уравнения второго порядка решение уравнения по формуле Эйлера в Microsoft Excel (таблички в режимах отображения чисел и формул графики точного и численного решений максимальную погрешность – абсолютную и относительную описание функций решения обыкновенного дифференциального уравнения odesolve и rkfixed; решения, полученные в пакете Mathcad с использованием функций odesolve и rkfixed. ВАРИАНТ УРАВНЕНИЕ НАЧАЛЬНЫЕ УСЛОВИЯ ПРОМЕЖУТОК ИНТЕГРИРОВАНИЯ Вариант 1. 0 6 5 4 y y y 4 , 0 ) 2 0 ( y ; 5 , 0 ) 2 0 ( y x [0.2, 2.5]. Вариант 2. 0 8 2 y y y 85 , 0 ) 4 0 ( y ; 5 , 1 ) 4 0 ( y x [0.4, 2.5]. Вариант 3. 0 7 8 y y y 5 , 0 ) 1 0 ( y ; 5 , 1 ) 1 0 ( y x [0.1, 2.5]. Вариант 4. 0 4 y y y 85 , 0 ) 3 0 ( y ; 5 , 1 ) 3 0 ( y x [0.3, 2.1]. Вариант 5 0 4 9 5 y y y 5 ) 1 0 ( y ; 8 2 ) 1 0 ( y x [0.1, 2.83]. Вариант 6. 0 34 2 y y 81 , 0 ) 3 0 ( y ; 5 , 2 ) 3 0 ( y x [0.3, 2.1]. Вариант 7. 0 8 4 y y 1 ) 8 0 ( y ; 5 , 0 ) 8 0 ( y x [0.8, 1.8]. Вариант 8. 0 2 y y 85 , 0 ) 1 ( y ; 5 , 1 ) 1 ( y x [-1, 0.65]. Вариант 9. 0 4 4 y y y 85 , 1 ) 5 0 ( y ; 5 , 0 ) 5 0 ( y x [-0.5, 1.8]. Вариант 10. 0 9 y y 8 , 1 ) 6 0 ( y ; 5 , 1 ) 6 0 ( y x [0.6, 1.9]. ТЕМА 8. ЧИСЛЕННОЕ РЕШЕНИЕ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ВТОРОГО ПОРЯДКА. КРАЕВАЯ ЗАДАЧА Задание 15. Решить первую краевую задачу для обыкновенного дифференциального уравнения x y y на промежутке изменения аргумента x [0, 1] шагом h=0,2. Краевые условия имеют вид y(0)=y(1)=0. Решение. В данной задаче x 0 =0; x 1 =0.2; x 2 =0.4; x 3 =0.6; x 4 =0.8; x 5 =1; p(x i )=0; q(x i )=-1; f(x i )=x i, y(x 0 )=0; y(x 5 )=0. Значит, система разностного метода будет иметь вид 0 04 , 0 ) ( 0 2 2 , 0 1 ) ( ) 1 ( 04 , 0 2 ) ( 0 2 2 , 0 1 04 , 0 ) ( 0 2 2 , 0 1 ) ( ) 1 ( 04 , 0 2 ) ( 0 2 2 , 0 1 04 , 0 ) ( 0 2 2 , 0 1 ) ( ) 1 ( 04 , 0 2 ) ( 0 2 2 , 0 1 04 , 0 ) ( 0 2 2 , 0 1 ) ( ) 1 ( 04 , 0 2 ) ( 0 2 2 , 0 1 0 ) ( 5 4 5 4 3 3 4 3 2 2 3 2 1 1 2 1 0 0 x y x x y x y x y x x y x y x y x x y x y x y x x y x y x y x y (8.1) Произведя вычисления в скобках перед неизвестными и подставив значения x i , можем записать систему в виде 04 , 0 8 , 0 04 2 04 , 0 6 , 0 04 2 04 , 0 4 , 0 04 2 04 , 0 2 , 0 04 2 4 3 4 3 2 3 2 1 2 1 x y x y x y x y x y x y x y x y x y x y (8.2) Эту систему четырех линейных алгебраических уравнений относительно четырех неизвестных y(x 1 ); y(x 2 ), y(x 3 ); y(x 4 ), можем решить в пакете Microsoft Excel методом итерации. Для этого выражаем в каждом уравнении неизвестные величины, стоящие на главной диагонали через остальные неизвестные 04 2 04 , 0 8 , 0 04 2 04 , 0 6 , 0 04 2 04 , 0 4 , 0 04 2 04 , 0 2 , 0 3 4 4 2 3 3 1 2 и выполняем вычисления до достижения требуемой точности. В первую строку таблицы Microsoft Excel записываем начальные приближения неизвестных (можем принять за нулевое приближение нулевые значения. В соответствующих столбцах следующей строки набираем формулы для вычисления переменных и копируем их до достижения решением требуемой точности (рис. 8.1). Рис. 8.1. Решение СЛАУ задачи 15 методом итерации На рис. 8.2 фрагмент таблицы с решением в режиме отображения формул. Рис. 8.2. Фрагмент таблицы с решением методом итерации в режиме отображения формул Для систем линейных алгебраических уравнений с трехдиагональными матрицами коэффициентов применим метод прогонки, являющийся модификацией метода Гаусса. Заполним таблицу для выполнения вычислений следующим образом. Впервой строке надпишем столбцы именами переменных, которые будут находиться в этих столбцах i – номер неизвестной f i -значения столбца свободных членов уравнений системы a i , b i , c i – значения коэффициентов в уравнениях системы (a i - коэффициенты перед неизвестными ниже главной диагонали, b i - коэффициенты перед неизвестными на главной диагонали, c i - коэффициенты перед неизвестными выше главной диагонали s i , g i , x i – столбцы, в которых будут проводиться вычисления. В столбцах F, G таблицы Microsoft Excel выполняем вычисления прямого хода. Как обычно формулы набираются во второй строке вычислений и далее копируются для всех уравнений системы. В столбце H выполняем вычисления обратного хода. Причем начинаем вычисления с последнего значения системы, далее пишем формулу для определения предпоследнего значения и копируем ее, захватив за правый нижний угол ячейки в ячейки предыдущих строк столбца таблицы. На рис. 8.3 приведено вычисление решения системы методом прогонки. Рис. 8.3. Решение методом прогонки в режиме отображения чисел Рис. 8.5. Решение методом прогонки в режиме отображения формул Полученное решение представим графически. Добавляем граничные условия и строим график (рис. 8.6). Рис. 8.6. График решения краевой задачи Можно ли проверить справедливость полученных результатов Расчетные значения определялись, как решения разностных уравнений. Следовательно, подставив в любое из уравнений системы соответствующие вычисленные значения y i , должны получить тождество. Такую проверку выполним для первого и второго неизвестных – данные в ячейках J2, J3 на рис. 8.7. Подставим в разностные уравнения полученные значения неизвестных. Как видим, они точно совпадают с правой частью этих уравнений, которые находятся в ячейках В, В. 1 2 3 4 5 A B C D E F G H i fi ai bi ci si gi xi 1 -0,0039 0 1 -0,49 -0,49 -0,0039 -0,0286 2 0,016 1 -2,04 1 -0,645 -0,0129 -0,0503 3 0,024 1 -2,04 1 -0,717 -0,0264 -0,0581 4 0,032 1 -2,04 0 -0,0442 -Решение краевой задачи первого рода 0 0 0,2 0,4 0,6 0,8 1 |