Ю11.3

Задача. Метод Эйткена-Стеффенсона (развитие метода простой итерации). Найти решение уравнения [latex]x=\varphi(x)[/latex] методом Эйткена-Стеффенсона, в котором от заданного начального [latex]x_{0}[/latex] три очередных приближения находятся по формулам: [latex]x_{n+1}=\varphi (x_{n})[/latex]; [latex]x_{n+2}=\varphi (x_{n+1})[/latex]; [latex]x_{n+3}=\frac{x_{n}x_{n+2}-x_{n+1}^{2}}{x_{n}-2x_{n+1}+x_{n+2}}[/latex]. Процесс продолжается до достижения одного из условий: [latex]\left | x_{n+3}-x_{n+2} \right |<\varepsilon[/latex] или [latex]x_{n}-2x_{n+1}+x_{n+2}=0[/latex].

Решение. В начале мы задаем начальное приближение  [latex]x=x_{0}[/latex] и [latex]\varepsilon[/latex]. Нам понадобится цикл с постусловием. Далее вычисляем первое приближение [latex]x_{1}=\varphi (x_{0})[/latex] и второе приближение [latex]x_{2}=\varphi (x_{1})[/latex]. Находим по формуле [latex]x=\frac{x_{0}x_{2}-x_{1}^{2}}{x_{0}-2x_{1}+x_{2}}[/latex]. Далее мы проверяем выполнение условия выхода из цикла: если [latex]\left | x-x_{0} \right |>\varepsilon[/latex] (нужно проверять разность [latex]\left | x-x_{0} \right |>\varepsilon[/latex], а не [latex]\left | x-x_{2} \right |>\varepsilon[/latex]) и [latex]x_{0}-2x_{1}+x_{2}\neq 0[/latex] тогда цикл продолжается, но если хотя бы одно из условий выполнено, то мы получаем, что [latex]x[/latex] — решение уравнения [latex]x=\varphi(x)[/latex].

Данная программа написана для функции [latex]9.2x^{3}+3.3x^{2}+4x-1[/latex].

Тесты:

x [latex]\varepsilon[/latex] Корень уравнения Комментарий
2 0.005 1.99889 Пройден
1 0.00002 0.233905 Пройден
-4 0.0005 -3.99981 Пройден

Код программы:

Код программы можно проверить здесь.

Решение на Java:

Ссылка на решение.

Ю11.8

Метод Симпсона.  Вычислить определённый интеграл [latex]I=\int_{a}^{b}{f\left(x \right)}dx[/latex] по формуле Симпсона: [latex]I\approx \frac{b-a}{6n}\left(y_{0}+4y_{1}+2y_{2}+\cdot \cdot \cdot +4y_{2n-1}+y_{2n} \right)[/latex], где [latex]2n[/latex] — количество отрезков разбиения, [latex]y_{0}[/latex], [latex]y_{1}[/latex], …, [latex]y_{2n}[/latex] — значение функции [latex]f\left(x \right)[/latex] на концах отрезков.
В задачах на численное интегрирование определённый интеграл требуется найти с заданной точностью, для чего вычисление по формуле метода рекомендуется проводить многократно, каждый раз уменьшая шаг интегрирования в два раза, пока разница между соседними приближениями не станет меньше заданной погрешности.

Функция [latex]a[/latex] [latex]b[/latex] [latex]eps[/latex] Интеграл Комментарий
[latex]f\left(x \right)=\sin \left(x^{2}+2x \right)[/latex] 1 3 0.0001 -0.143058 Тест пройден.
[latex]\ln \left(1+x \right)[/latex] 1 3 0.0001 2.15888 Тест пройден.
[latex]\tan \left(3x^{3} \right)[/latex] 2 15 0.01 0.0256033 Тест пройден.
[latex]x\left( x^{2}-1\right)\left(x+1 \right)[/latex] -1 1 0.3 -0.265625 Тест пройден.

C++:

Java:

С помощью программы можно вычислить интеграл любой непрерывной на  [latex]\left[a;b \right][/latex]  функции, для этого нужно изменить 7 строку.

В условии главного цикла [latex]N\leq 4[/latex], так как важно, чтобы цикл выполнился хотя бы дважды. Потому что первый раз мы сравниваем новое значение интеграла не с предыдущим вычисленным, а с нулём, и если новое значение интеграла будет меньше погрешности, то цикл прекратится после первого же выполнения (без условия [latex]N\leq 4[/latex]).

Задача на Ideone:
C++
Java

Ю11.13

Задача

Метод Рунге-Кутта. Найти приближенное решение обыкновенного дифференциального уравнения  [latex]y^\prime=f(x,y), y(a)=y_{0}[/latex] методом Рунге-Кутта пятого порядка на отрезке [latex][a,b][/latex] с заданным шагом [latex]h[/latex]. Значения функции [latex]y(x)[/latex] в узловых точках вычисляется по формуле: [latex]y_{i+1}=y_{i}+\frac{h}{6}(k_{1}+2k_{2}+2k_{3}+k_{4}), i=0,1,2,\cdots[/latex], где [latex]k_{1}=f(x_{i},y_{i}); k_{2}=f(x_{i}+\frac{h}{2},y_{i}+\frac{h}{2}k_{1});[/latex][latex]k_{3}=f(x_{i}+\frac{h}{2},y_{i}+\frac{h}{2}k_{2}); k_{4}=f(x_{i}+h,y_{i}+hk_{3})[/latex].

Решим дифференциальное уравнение такого вида: [latex]y^\prime=x+y[/latex] при начальном условии [latex]y(0)=1[/latex] на отрезке [latex][0, 0.5][/latex] с шагом интегрирования [latex]h=0.1[/latex]

Screenshot_1

Screenshot_2

 

Screenshot_3

Код программы

Для решения примера введем данные во входной поток в таком порядке: [latex]0.0[/latex], [latex]0.5[/latex], [latex]0.1[/latex], [latex]1.0[/latex], где первое и второе число — начало и конец отрезка интегрирования соответственно, третье — шаг интегрирования и четвертое — значение [latex]y[/latex] в точке [latex]a[/latex] (в начале отрезка), т.е. [latex]y(a)=y(0)[/latex].

В программе присутствует функция, которой мы передаем параметры [latex]x, y[/latex] и которая возвращает само дифференциальное уравнение. Далее, в цикле высчитываем значения [latex]k_{1},k_{2},k_{3},k_{4}[/latex], передавая каждый раз параметры в функцию с шагом [latex]h[/latex] до тех пор пока не дойдем до конца промежутка. После завершения цикла выводим значение [latex]y_{0}[/latex].

 

Ю11.6

Задача. Метод прямоугольников. Вычислить определенный интеграл [latex]I=\int_{a}^{b}{f(x)dx}[/latex] методом прямоугольников: [latex]\int^b_a f(x)\,dx \approx h (\frac{y_0}{2} + y_1 + \ldots + y_{n-1}+\frac{y_{n}}{2})[/latex], где [latex]n[/latex] — количество отрезков разбиения;  [latex]y_{0},y_{1},…,y_{n}[/latex] — значения функции на концах отрезков.

Вычислим для функции [latex] f(x)=2x^{3}-7x+4[/latex]:

[latex] \int_{0}^{2}{(2x^{3}-7x+4)dx}=2[/latex]

Решение: 

Введена функция, которая подсчитывает значение в точке.  Согласно формуле в условии, вычисляем требуемое значение.

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

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

Количество отрезков разбиения на [a,b] Результат
50 2.0032
500 2.000032
1000 2.000008
5000 2.00000032

С работой программы можно ознакомиться здесь.

Ю11.7

Метод трапеций. Вычислить определенный интеграл [latex]I=\int_{b}^{a}f(x)dx [/latex] методом трапеций:[latex] I\approx \frac{b-a}{2n}(y_{0}+2y_{1}+\dots+2y_{n-1}+y_{n}), [/latex] где [latex] n [/latex] — количество отрезков разбиения; [latex]y_{0},y_{1},\ldots,y_{n} [/latex] — значения функции [latex]f(x) [/latex] на концах отрезков.

Вычислим определенный интеграл для функции [latex]y=-3x^2+2x+9[/latex]
[latex] \int_{-1}^{2}(-3x^2+2x+9)dx=21 [/latex]

Решение:

Ссылка на ideone C++: http://ideone.com/RJpYSw

Ссылка на ideone Java: http://ideone.com/AfEDeq

 

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

Ю11.15

Метод парабол. Найти минимум заданной функции [latex]y=f(x)[/latex], двигаясь от заданной точки [latex]x_{0}[/latex] по методу парабол:

[latex]x_{i+1}=x_{i}-\frac{h}{2}\frac{f\left( x_{i}+h\right)-f\left( x_{i}-h\right)}{f\left( x_{i}+h\right)-2f\left( x_{i}\right)+f\left( x_{i}-h\right)}, i=0,1,\ldots[/latex], пока не будет достигнута заданная точность.

Функция [latex]x^{3}+10\sin (5x)[/latex]:

Безымянный

[latex]x_{0}[/latex] [latex]\varepsilon[/latex] Точка минимума по графику. Точка минимума по программе. Комментарий.
0 0,001 -0,315353 -0,315353 Тест пройден.
1 0,0001  0,932048 0,932048 Тест пройден.
2 0,0001  2,14327 2,14327 Тест пройден.
-3 0,01 -2,93616 -2,93616 Тест пройден.
-2 0,0001  -1,6017 -1,6017 Тест пройден.

Код программы (C++):

Для функции [latex]x^{3}+10\sin (5x)[/latex] (функцию всегда можно изменить, достаточно исправить строку 6):

Java:

 

Чтобы сделать программу, я почитала “Численные методы” Н. Н. Калиткина, где и было сказано, что в качестве вспомогательного шага [latex]h[/latex] при расчётах на ЭВМ обычно выбирают значение 0,001.

Далее пользователю предоставляется возможность ввести значение [latex]x_{0}[/latex] и задать точность [latex]\varepsilon[/latex].  Корень [latex]x[/latex] является минимумом функции тогда и только тогда, когда вторая производная этой функции больше 0. Поэтому мы запускаем цикл, который будет сдвигать наше [latex]x_{0}[/latex] в положительном направлении оси [latex]x[/latex], пока вторая производная функции не будет удовлетворять условию задачи.

Далее мы вычисляем значение [latex]x_{i+1}[/latex], и запускаем цикл, который будет продолжаться до тех пор, пока разница между [latex]x_{i}[/latex] и [latex]x_{i+1}[/latex]  не будет меньше заданной погрешности [latex]\varepsilon[/latex].

Затем на экран выводится сама точка минимума функции.

Программу можно посмотреть здесь (C++) и здесь (Java).

Ю3.33

Задача: В задаче задана функция и её разложение в ряд или произведение. Численно убедится в справедливости равенства, для чего для заданного значения аргумента [latex] x [/latex] вычислить левую его часть и разложение, стоящее в правой части, с заданной погрешностью [latex]\varepsilon\[/latex]. Испытать разложение на сходимость при разных значениях аргумента, оценить скорость сходимости, для чего вывести число итераций [latex] n [/latex]  (слагаемых или сомножателей), необходимых для достижения заданной точности.
[latex]\sin \left(x \right) = x \left(1 — \frac{x^{2}}{\pi^{2}} \right) \left( 1 — \frac{x^{2}}{4\pi^{2}} \right) \ldots \left(1 — \frac{x^{2}}{ \left( n-1 \right)^{2}\pi^{2}} \right)[/latex]

Тесты:

x n-Excel sin — Excel deviation — Excel e — input exact sinus n — program deviation — program
12  22  -1,029367  0.56  0.56  -0.536573  22  0.524789
5  7  -1,346966
 0.54  0.54  -0.958924  7  0.461469
2.5  2  0,771704  -1.15  -1.15  0.598472  2  -0.318384
1  2  0,875915  -1.38  -1.38  0.841471  2  -0.057208
0  2  0  -0.53  -0.53  0  2  0
-1  2  -0,875915  0.3  0.3  -0.841471  2  0.057208
-2.5  6  -0,659746  0.08  0.08  -0.598472  6  0.073087
-5  4  1,700161  -1.62  -1.62  0.958924  4  -1.061024
-12  12  1,755690  -1.63  -1.63  0.536573  12  -1.417062

Код программы:

Код программы на языке Java:

Ссылка:http://ideone.com/nZluY8

Программа состоит из следующих частей:

  1. Определение вспомогательной функции
  2. Рабочие переменные
  3. Определяем лимит числа итераций чтобы избежать бесконечного цикла
  4. Ввод данных — аргумента функции и заданной погрешности
  5. Условие выполнения
  6. Вывод результата

Правильность результата проверялась на калькуляторе. тестирование показало, что правильные значения получаются только при очень большом  n  .  Только начиная с десятков тысяч результат близок к показанию калькулятора.

Данные для проверки подготавливались в Excel ввиду относительно большого количества итераций, необходимых для расчётов. Полученная погрешность с точным значением синуса использовались как входное [latex]e[/latex]. Результирующее количество итераций программы сравнивалось с числом итераций, сделанных в Excel. Ввиду того, что double намного точнее чем значения в Excel, возможны различия в числе итераций ( назначенных ).

Ссылка на ideone.com: http://ideone.com/fork/Wh91nH

 

Ю11.16

Задача:
Для заданной матрицы [latex]A(n, n)[/latex] найти обратную [latex]A^{-1}[/latex], используя итерационную формулу:
[latex]A_{k}^{-1} = A_{k-1}^{-1}(2E-A A_{k}^{-1}),[/latex]
где [latex]E[/latex] — единичная матрица, [latex]A_{0}^{-1} = E[/latex]. Итерационный процесс заканчивается, если для заданной погрешности [latex]\varepsilon[/latex] справедливо:
[latex]\left| det(A A_{k}^{-1})-1 \right| \le \varepsilon[/latex]

Анализ задачи:
Прежде чем приступать к решению средствами языка C++, я создал прототип в системе компьютерной алгебры Wolfram Mathematica, с которым планировал сверяться при тестировании программы. Тем не менее, внимательный анализ показал, что с таким выбором начального приближения процесс уже на пятом шаге в значительной мере расходится даже для матриц размера 2*2. После уточнения условий и анализа дополнительного материала, посвященного методу Ньютона-Шульца, исходное приближение было изменено (по результатам исследования «An Improved Newton Iteration for the Generalized Inverse of a Matrix, with Applications», Victor Pan, Robert Schreiber, стр. 8):
[latex]A_{0}^{-1} =\frac { A }{ \left\| A \right\|_{1} \left\| A \right\| _{\infty } }[/latex], где [latex]{ \left\| A \right\| }_{ 1 }=\max _{ i }{ \sum _{ j=0 }^{ n-1 }{ \left| { a }_{ ij } \right| } } [/latex], [latex]{ \left\| A \right\| }_{ \infty }=\max _{ j }{ \sum _{ i=0 }^{ n-1 }{ \left| { a }_{ ij } \right| } }[/latex].
Эффективность предложенного подхода иллюстрируют результаты работы прототипа:
процесс сходится
Следовательно, из пространства задачи можно переместиться в пространство решения и составить алгоритм реализации предложенного метода на языке C++.

Тесты:

[latex]n[/latex] [latex]A[/latex] [latex]A^{-1}[/latex] Результат
3 1 2 3
5 5 7
11 13 7
-0.964771 0.430661 -0.017183
0.723533 -0.447973 0.137884
0.172358 0.155200 -0.086211
Пройден
2 1 2
2 1
-0.333067 0.666400
0.666400 -0.333067
Пройден
5 1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
Пройден
3 1 2 3
4 5 6
7 8 9
Матрица вырождена Пройден

Программный код:

Программа доступна для тестирования по ссылке: http://ideone.com/7YphoX.

Алгоритм решения
При решении данной задачи оправдывает себя подход «разделяй и властвуй»: вычисление обратной матрицы подразумевает промежуточные этапы, корректной реализации которых следует уделить особое внимание. Наверняка стоило бы написать класс matrix и реализовать перегрузку операций, но задачу можно решить и без применения объектно-ориентированных средств, пусть от этого решение и потеряет в изящности.
Можно выделить следующие подзадачи:

  1. Инициализация динамического двухмерного массива и передача его в функцию.
  2. Вычисление определителя матрицы (с применением метода Гаусса).
  3. Вычисление нормы матрицы.
  4. Транспонирование матрицы.
  5. Умножение матрицы на скаляр.
  6. Матричное умножение матриц.
  7. Сложение двух матриц.
  8. Непосредственно итерационный процесс Ньютона-Шульца

Ниже приведено пояснение к подходу к реализации некоторых подзадач:

  • Выделение памяти для хранения массива происходит не на этапе запуска программы, после компиляции. Для инициализации использован конструктор new
  • Вычисление определителя можно разбить на два последовательных шага:
    1. Приведение матрицы к верхнетреугольному виду (методом Гаусса).
    2. Вычисление определителя как произведения элементов главной диагонали.
    3. Если матрица вырождена, то дальнейшие вычисления не производятся.

  • Нормы матрицы вычисляются как максимальные значения суммы элементов по столбцам и строкам соответственно.
  • При транспонировании обмениваются местами элементы, симметричные главной диагонали.
  • При умножении матрицы на скаляр каждый элемент матрицы умножается на соответствующее вещественное число.
  • При перемножении двух квадратных матриц используется промежуточный массив для хранения результата вычислений.
  • Сложение двух матриц аналогично попарному сложению элементов, расположенных на соответствующих позициях.
  • Максимально допустимая погрешность для метода Ньютона-Шульца [latex]\varepsilon = 0.001[/latex]. Программа использует локальный массив для хранения [latex]A_{k}^{-1}[/latex], инициализация которого происходит в теле цикла.

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

Ю11.12

Задача:
Интерполяционный многочлен Лагранжа. Значения функции [latex]y=f\left(x\right)[/latex] заданы таблично в массиве [latex]Y\left(x\right)[/latex] при соответствующих значениях аргумента в упорядоченном массиве [latex]X\left(x\right)[/latex]. Найти значение функции в произвольной точке [latex]x[/latex] по формуле Лагранжа:
[latex]y={L}_{n}\left(x\right)=\sum _{i=1}^{n}{{y}_{i}\prod _{\underset{j\neq i}{j=1}}^{n}{\frac{x-{x}_{j}}{{x}_{i}-{x}_{j}}}}[/latex]

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

X: -1 0 1 2
Y: 1 2 3 -1

Код программы:

Код программы на Java:

 

Идея решения:
Использование отдельной функции, которой надо передать 3 параметра: значение аргумента, массив исходных значений аргумента и массив исходных значений функции при соответствующих аргументах. Эта функция возвращает значение многочлена Лагранжа при подстановки вместо [latex]x[/latex] конкретного значения.

В функции был использован внешний цикл for для вычисления суммы в формуле и внутренний цикл for для вычисления произведения [latex]\prod _{\underset{j\neq i}{j=1}}^{n}{\frac{x-{x}_{j}}{{x}_{i}-{x}_{j}}}[/latex].

В программе мы разбиваем отрезок [-3, 3] на 101 отдельную точку и вычисляем значение полинома Лагранжа для каждой из этих точек. Выводим на экран.

График:ГрафикПолиномаЛагранжа

Комментарии: Точки выбраны специально. Мне просто захотелось увидеть, как выглядит график функции, которую мы разбирали на паре (кстати по той же самой теме =) ). График сделан в xls. Векторы выбраны по причине ограниченности массива.

А58г

 Задача: Дано действительное число  [latex]a[/latex]. Для функции  [latex]f(x)[/latex], график которой представлен на рисунке, вычислить  [latex]f(a)[/latex].
График:
a
Тесты:

a f(a)
1 1
3.2 -0.015371
6 -0.027469
0 0
-1 1
-2.5 2.5
1.5 1
1.8 1
1.001 1

Код программы:

Код программы на языке Java:

Ссылка:http://ideone.com/e6UFys

Результат вычисляем по формуле:
[latex]y = ka + b[/latex]
Программа состоит из следующих частей:

  1. Объявление переменных a, y, k, b типа float для хранения данных
  2. Ввод пользователем значений переменной а с помощью scanf
  3. Вычисление и вывод результата по формуле с предварительным сравнением значения а
  4. Завершение программы

Программа сравнивает значение переменной [latex]a[/latex] с значениями переменной [latex]x[/latex] на четырёх диапазонах, и в зависимости от диапазона использует для функции [latex]y = ka + b[/latex] нужные значения [latex]k[/latex] и [latex]b[/latex]. Так вычисляется [latex]f(a)[/latex].
Ссылка на ideone.com : http://ideone.com/N2toyp

Ю2.10

Задача: Посылка. Можно ли коробку размером  [latex]a \times b \times c[/latex]  упаковать в посылку размером  [latex]r \times s \times t[/latex] ? «Углом» укладывать нельзя.

Тесты:

a b c r s t result
1 1 1 3 3 3 Package is fitting to the box
23 2 4 54 45 22 Package is fitting to the box
0 0 0 1 1 0 Package cannot fit !
109 122 222 11 22 33 Package cannot fit !
13 43 21 55 76 89 Package is fitting to the box
23 15 17 44 81 92 Package is fitting to the box

 

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

 

Код программы на языке Java:

 

Ссылка http://ideone.com/cRnZc7

Программа состоит из следующих частей:

  1. Объявление переменных a,b,c,r,s,t типа float для хранения входных данных
  2. Ввод пользователем значений переменных a, b, c, r, s, t с помощью scanf
  3. Вывод исходных данных
  4. Сравниваем рёбра пакета с каждым ребром коробки
  5. Завершение программы

 

Программа сверяет размеры коробки ( package ) с размерами посылки ( box ). Если размеры коробки меньше, чем посылки, программа сообщает о возможности упаковки коробки в посылку. В противном случае выводиться сообщение о невозможности упаковки коробки данных размеров в посылку.

Ссылка на ideone.com :  http://ideone.com/XYcV5W