Пример курсовой работы по прикладной математике: Вычислительная математика

3. Аппроксимация функций

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

Одним из распространенных и практически важных случаев связи между аргументом и функцией является задание этой связи в виде некоторой таблицы {x_i ; y_i}, например, экспериментальные данные. На практике часто приходится использовать табличные данные для приближенного вычисления у при любом значении аргумента х (из некоторой области). Этой цели служит задача о приближении (аппроксимации) функций: данную функцию f(x) требуется приближенно заменить некоторой функцией g(х) так, чтобы отклонение g(х) от f(x) в заданной области было наименьшим. Функция g(х) при этом называется аппроксимирующей. Если приближение строится на заданном дискретном множестве точек {x_i}, то аппроксимация называется точечной. К ней относятся интерполирование, среднеквадратичное приближение и др. При построении приближения на непрерывном множестве точек (например, на отрезке [a, b]) аппроксимация называется непрерывной или интегральной. MATLAB имеет мощные средства точечной и непрерывной аппроксимации с визуализацией результата.

Рисунок 2-алгоритм формирования матрицы Грама и аппроксимации полиномом

3.2 Примеры

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

Таблица 1

Оценить погрешность эмпирических формул.

Протокол программы. В окне команд набираются значения x_i и y_i. Далее выполняется команда построения графика только узловых точек.

>>x = [0.75, 1.50, 2.25, 3.00, 3.75] ;

>>y = [2.50, 1.20, 1.12, 2.25, 4.28] ;

>>рlot(x, y, 0 ) ;

Появляется окно с символами 0 на месте узловых точек.

Рисунок 3

Внимание. Следует помнить, что при полиномиальной аппроксимации максимальная степень полиноме на 1 меньше числа экспериментальных точек!

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

— показать уравнение аппроксимирующей функции у = g(х);

— выбрать метод подбора: сплайн интерполяции

эрмитовская интерполяция

линейный

квадратные

кубический

…………….

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

Рисунок 4

Чтобы узнать погрешность аппроксимации, надо отметить птичкой параметр График остатка в окне Основной Монтаж, и Показать норму остатков. График погрешностей с нормами можно вынести в отдельное окно, или вместе с графиком и аппроксимирующих функций — суб-график. Норма погрешностей указывает на статистическую оценку среднеквадратической погрешности. Чем она меньше, тем точнее полученная аппроксимирующая функция у = g(х). В нашем примере:

Linear : norm of residuals (норма погрешности) = 2.1061

Quadratic : norm of residuals = 0.10736

Cubic : norm of residuals = 0.035857

4 th degree : norm of residuals = 9.6305e-015.

График погрешностей можно выводить в виде диаграмм (зоны), линий (линии) или отдельных точек (фрагменты). Сам график погрешностей представляет собой зависимость разности g(х) f(x) в условных точках, соединенных прямыми линиями.

Кроме линейной и полиномиальной аппроксимации можно выбрать сплайн аппроксимацию — когда на каждом интервале приближения используется кубический полином с новыми коэффициентами. В этом случае нельзя получить выражение для аппроксимирующей функции, т. е. такая аппроксимация является неполной. Аналогичными свойствами обладает и Эрмитовая аппроксимация. Она имеет только графическую интерпретацию.

4. Численное решение обыкновенных дифференциальных уравнений

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

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

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

Требуется найти функцию у = у(х), удовлетворяющую уравнению

алгебраический аппроксимация дифференциальный интеграл

и принимающую при х = х₀ заданное значение у₀:

. (2)

При этом решение необходимо получить в интервале х₀хх_к. Из теории дифференциальных уравнений известно, что решение у(х) задачи Коши (1), (2) существует, единственно и является гладкой функцией, если правая часть F(x, y) удовлетворяет некоторым условиям гладкости. Численное решение задачи Коши методом Рунге-Кутта 4-го порядка заключается в следующем. На заданном интервале [х₀, х_к] выбираются узловые точки. Значение решения в нулевой точке известно у(х₀) = у₀. В следующей точке у(х₁) определяется по формуле

, (3)

Здесь

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

Рисунок 5 — схема метода Эйлера-Коши с уточнением

4.2 Примеры

Решить задачу Коши

y’ = f ,

f’ = (3xі + x·f)/ хІ,

u’ = v,

v’ = v/ х,

y(1) = 2, f(1) = б, u(1) = 0, v(1) = 1.

. (1)

Точное решение имеет вид

Выполним решение данной задачи с помощью программы ode45. Вначале в М-файл записываем правую часть уравнения (1), сам М-файл оформляется как файл — функция, даем ему имя F:

function dydx = F(x, y)

dydx = zeros(1,1);

dydx(1) = 2*(x^2+y(1));

Для численного решения задачи Коши в окне команд набираются следующие операторы.

Протокол программы.

>>[XY] = ode45 ( @ F , [0 1] , [1] ) ;

% Дескриптор @ обеспечивает связь с файлом — функцией правой части

% [0 1] — интервал на котором необходимо получить решение

% [1] — начальное значение решения

>> рlot (X,Y) ;

>> % Построение графика численного решения задачи Коши (1)

>> hold on; gtext ( y(x) )

% Команда позволяет с помощью мышки нанести на график надпись у(х)

>> [XY]

>> % Последняя команда выводит таблицу численного решения задачи.

Результаты решения. График решения задачи Коши (1) показан на рис. 2. Численное решение представлено в таблице 2, где приведены только отдельные узловые точки. В программе ode45 по умолчанию интервал разбивается на 40 точек с шагом h = 1/40 = 0.025.

Рисунок 6

Таблица 3

Как следует из таблицы 3 численное решение программой ode45 является точным.

5. Приближенное вычисление определенных интегралов

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

К вычислениям определенных интегралов сводятся многие практические задачи физики, химии, экологии, механики и других естественных наук. На практике формулой Ньютона-Лейбница не всегда удается воспользоваться. В этом случае используются методы численного интегрирования. Они основаны на следующих соображениях: с геометрической точки зрения определенный интеграл представляет собой площадь криволинейной трапеции. Идея методов численного интегрирования сводится к разбиению интервала [a; b] на множество меньших интервалов и нахождению искомой площади как совокупности элементарных площадей, полученных на каждом частичном промежутке разбиения. В зависимости от использованной аппроксимации получаются различные формулы численного интегрирования, имеющие различную точность. Рассмотрим методы трапеций и Симпсона (парабол).

5.2 Метод трапеций

Здесь используется линейная аппроксимация, т. е. график функции y = f(x) представляется в виде ломаной, соединяющей точки y_i.

Формула трапеций при постоянном шаге , где п — число участков, имеет вид

. (2.8)

В MATLAB данную формулу реализуют программы

cumtrapz(у)

cumtrapz(x,y)

trapz(y)

trapz(x,y)

Наиболее интересна последняя из них. Данная функция вычисляет интеграл от функции у(х) по х методом трапеций. Аргумент и функция задаются в виде векторов или х — в виде вектора, а у — в виде матрицы. Если у(х) — матрица, то функция возвращает вектор значений интеграла для каждого столбца матрицы.

5.3 Метод Симпсона

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

. (2.9)

В MATLAB формула Симпсона реализуется программой quad(name, a, b [,tol, trace]), где

name — имя М-функции, задающей подынтегральное выражение;

a, b — пределы интегрирования;

tol — точность вычисления;

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

Подынтегральная функция может задаваться с помощью дескриптора @, тогда она программируется в файле — функции, или с помощью апострофов, тогда она записывается в самой программе quad. Точность вычисления интегралов по умолчанию принята равной 110^-6.

Рисунок 7 — алгоритма метода Симпсона

5.4 Вычисление интегралов

Символьные вычисления неопределенных интегралов в MATLAB осуществляется при помощи функции:int(fun, var),где fun — символьное выражение, представляющее собой подынтегральную функцию, а var — переменная интегрирования.

Пример вычисления неопределенного интеграла:

symsx %Определение символьной переменной

f=sym(‘exp(x) -x’); %Определение символьной функции

int(f,x) %Вычисление неопределенного интеграла

Результатом будет:

ans =

exp(x)-1/2*x^2

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

int(fun, var, a, b),

где fun -подынтегральная функция, а var — переменная интегрирования, a, b — пределы интегрирования.

Пример вычисления определенного интеграла:

I1=int(‘exp(x)-x’,’x’,-1,0); %Символьное решение

vpa(I1,5)%Численное решение

Результатом будет:

ans =

1.1321

5.5 Примеры

Пример: Вычислить и вывести на печать по методам трапеций и Симпсона значения интеграла

Протокол программы метода трапеций

>>x = 0 : 0.0001 : 1.0 ;

>>y = 1./ (1+x.^2) ;

>>z = trapz(x, y)

Результат вычислений

0.7854

Протокол программы метода Симпсона

>>quad( (1./(1+x.^2)) , 0, 1)

Результат вычислений

ans = 0.7854

Точное значение интеграла равно 0.785398163.

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

6. Численное решение систем нелинейных уравнений

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

Задача нахождения корней нелинейных уравнений встречается в различных областях научно-технических исследований. Проблема формулируется следующим образом. Пусть задана непрерывная функция f(x) и требуется найти корень уравнения

f(x) = 0.

Будем предполагать, что имеется интервал изменения х [a; b], на котором необходимо исследовать функцию f(x) и найти значение х₀, при котором f(x₀) равно или весьма мало отличается от нуля.

Данная задача в системе MATLAB может быть решена следующим образом. Вначале необходимо построить график функции f(x) на заданном интервале и убедиться в существовании корня или нескольких корней. Затем применить программы поиска корней. Если существует один корень и график f(x) пересекает ось ох, то можно применить программу fzero. Если f(x) имеет больше одного корня и может касаться и пересекать ось ох, то следует применить более мощную программу fsolve из пакета OptimizationToolbox, которая решает задачу методом наименьших квадратов. Программа fzero использует известные численные методы: деление отрезка пополам, секущей и обратной квадратичной интерполяции.

Рисунок 8 — метод Ньютона для решения систем нелинейных уравнений

Так как метод Ньютона отличается высокой скоростью сходимости при выполнении условий сходимости, на практике критерием работоспособности метода является число итераций: если оно оказывается большим (для большинства задач >100), то начальное приближение выбрано плохо.

6.2 Примеры

Пример: Найти корень нелинейного уравнения 10^х + 2х — 100 = 0 на интервале [1.0; 2.0].

Протокол программы

>> % Строим график заданной функции

>>x = 1.0 : 0.001 : 2.0; y = 10.0.^x + 2.0*x — 100.0;

>>рlot (x, y) ; grid on

Появляется окно с графиком функции 10^х + 2х — 100 (см. рис. 1), из которого следует, что корень функции на заданном интервале существует. Для точного определения корня применяем fzeroиfsolve.

Рисунок 9

>>X1 = fzero ( (10.0.^x + 2.0*x — 100.0) , [1.0 2.0])

Результат решения

Х1 = 1.9824

>>X2 =fsolve ( (10.0.^x + 2.0*x — 100.0) , 1.0 : 2.0)

Результат решения

Х2 = 1.9824

7. Численное решение оптимизационных задач

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

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

, (1)

где х₁, х₂, … , х_п — проектные параметры.

Можно выделить 2 типа задач оптимизации — безусловные и условные. Безусловная задача оптимизации состоит в отыскании максимума или минимума функции (1) от п действительных переменных и определении соответствующих значений аргументов на некотором множестве Gn-мерного пространства. Обычно рассматриваются задачи минимизации; к ним легко сводятся и задачи на поиск максимума путем замены знака целевой функции на противоположный. Условные задачи оптимизации — это такие, при формулировке которых задаются некоторые условия (ограничения) на множестве G. Здесь рассмотрим только безусловные задачи оптимизации.

Поиск минимума функции одной переменной

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

Поиск минимума функций нескольких переменных

Данная задача значительно сложнее первой. Рассмотрим ее решение на примере функции двух переменных. Алгоритм может быть распространен на функции большего числа переменных. Для минимизации функций нескольких переменных MATLAB использует симплекс — метод Нелдера-Мида. Данный метод является одним из лучших методов поиска минимума функций многих переменных, где не вычисляются производные или градиент функции. Он сводится к построению симплекса в п-мерном пространстве, заданного п + 1 вершиной. В двухмерном пространстве симплекс является треугольником, а в трехмерном — пирамидой. На каждом шаге итераций выбирается новая точка решения внутри или вблизи симплекса. Она сравнивается с одной из вершин симплекса. Ближайшая к этой точке вершина симплекса заменяется этой точкой. Таким образом, симплекс перестраивается и позволяет найти новое, более точное положение точки решения. Алгоритм поиска повторяется, пока размеры симплекса по всем переменным не станет меньше заданной погрешности решения. Программу, реализующую симплекс-методы Нелдера-Мида, удобно использовать в следующей записи

[x, min f] = f min search ( … ),

где х — вектор координат локального минимума; minf — значение целевой функции в точке минимума. Саму целевую функцию удобно представить с помощью дескриптора @ в М-файле.

7.2 Примеры

Пример: Найти и вывести на печать минимальное значение функции (одной переменной)

f(x) = 24 — 2x /3 + x²/ 30 на [5; 20].

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

Протокол программы

>>x = 5.0 : 0.001 : 20.0 ; y = 24 — 2* x/3 + x.^2/30 ;

>>рlot(x, y) ; gridon

Появляется окно с графиком этой функции (рис. 8), где отмечаем наличие минимума.

Рисунок 10

Далее, для точного определения координаты и значения минимума привлекаем программу fminbnd.

>> [x, y] = fminbnd ( (24.0 — 2* x/3 + x.^2/30) , 5.0, 20.0)

Результат поиска

х = 10.0000

у = 20.6667

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

1. Мартынов Н.Н. Введение в MATLAB 6. -М.:КУДИЦ-ОБРАЗ, 2002.-352 с.

2. Мэтьюс Джон Г., Финк Куртис Д. Численные методы. Использование MATLAB. -М. Вильямс, 2001. — 720 с.

3. Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. М.: Наука, 1986. — 288 с.

4. Потемкин В.Г. Система инженерных и научных расчетов MATLAB 5.х: — в 2-х т. Том 1. — М.:ДИАЛОГ — МИФИ,1999. — 306 с.

5. Потемкин В.Г. Система инженерных и научных расчетов MATLAB 5.х: — в 2-х т. Том 2. — М.:ДИАЛОГ — МИФИ,1999. — 304 с.

6. Справочник по специальным функциям с формулами, графиками и таблицами / Под ред. М. Абрамовича и И. Стиган: Пер. с англ. М.: Наука, 1979. — 832 с.

7. Демидович Б.П., Марон И.А., Шувалова Э.З. Численные методы анализа. — М.:Наука, 1967. — 386 с.

8. Каханер Д., Моулер К., Нэш С. Численные методы и математическое обеспечение. — М.: Мир, 1998. — 575 с.

9. Копченова М.В., Марон И.А. Вычислительная математика в примерах и задачах. — М.: Наука, 1977. — 367 с.

Размещено на