Решение уравнений и интерполяция функций

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

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

Решение уравнений и интерполяция функций



Введение

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

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

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

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

1) описание математической модели задачи на основе физической или экономической модели;

2) изучение методов решения поставленной математической модели задачи и создание новых методов;

3) выбор метода решения задачи исходя из заданной точности решения и особенностей задачи;

4) составление блок-схемы программы для решения задачи на ЭВМ;

5) отладка программы и оценка полученных результатов;

6) решение задачи на ЭВМ, построение графиков, получение оценки погрешностей, обоснование результатов.

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



1. Решение нелинейных уравнений

Нелинейными уравнениями называются уравнения вида

(1)

где - нелинейная функция, которая может относиться к трем типам:

1) нелинейная алгебраическая функция вида

;

2) трансцендентные функции - тригонометрические, обратные тригонометрические, логарифмические, показательные и гиперболические функции;

3) различные комбинации этих функций, например.

Решением нелинейного уравнения (1) является такая точка, которая при подстановке в уравнение (1) обращает его в тождество. На практике не всегда удается подобрать такое решение. В этом случае решение уравнения (1) находят с применением приближенных (численных) методов. Тогда решением будет являться такая точка, при подстановке которой в уравнение (1) последнее будет выполняться с определенной степенью точности, т.е., где e - малая величина. Нахождение таких решений и составляет основу численных методов и вычислительной математики.

Решение нелинейных уравнений разделяется на два этапа: отделение корней уравнений и уточнение корней нелинейных уравнений.

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

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

Второй способ отделения корней нелинейных уравнений - аналитический. Он основывается на следующих трех теоремах.

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

Теорема 2. Если функция непрерывна на отрезке, выполняется условие вида и производная сохраняет знак на, то на отрезке имеется единственный корень.

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

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

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



1.1 Метод простых итераций

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

(2)

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

,

и учитывая непрерывность функции на отрезке, получим

или .

Корень можно вычислить с заданной точностью по итерационной формуле

(3)



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

(4)

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

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

а б

Рисунок 1.1. Сходящийся (а) и расходящийся (б) итерационные процессы

Часто, если итерационный процесс расходится из-за невыполнения условия (4), нелинейное уравнение можно привести к виду, допускающему сходящиеся итерации.

Выполнения условия сходимости можно добиться путем перехода от исходного уравнения к эквивалентному виду следующим образом: сначала умножить обе части уравнения (1) на, а затем прибавить к обеим частям x, тогда. Обозначив, получим уравнение (2). Константа с выбирается так, чтобы выполнялось достаточное условие сходимости итерационного процесса (4), т.е.

(5)

Условие (5) равносильно двойному неравенству

поэтому константа выбирается из соотношений:

Метод простых итераций и почти все другие итерационные методы имеют два достоинства:

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

2) позволяют достигнуть любой заданной точности при любом начальном приближении.

Недостатки методов:

- трудность приведения уравнения (1) к виду (2).

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

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

1) Когда два последних приближения отличается между собой по модулю на заданную величину e:

(6)

Одного критерия (6) недостаточно, так как в случае крутизны графика, данное условие будет выполнено, но может находиться далеко от корня;

2) Когда последнее вычисленное приближение к корню удовлетворяет уравнению с заданной точностью:

(7)

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

1.2 Метод Ньютона

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

.



Тогда итерационный процесс осуществляется по формуле:

(8)

Геометрически этот процесс представлен на рисунке 2. Он означает замену на каждой итерации k графика кривой касательной к ней в точках с координатами.

Рисунок 1.2. Геометрическая интерпретация метода Ньютона

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

(9)



При произвольном начальном приближении итерации сходятся, если

(10)

Метод Ньютона рекомендуется применять для нахождения простых действительных корней уравнения (1).

Достоинством метода является то, что он обладает скоростью сходимости, близкой к квадратичной.

Недостатки метода:

- не при любом начальном приближении метод Ньютона сходится, а лишь при таком, для которого;

- если, т.е. касательная к графику почти параллельна оси абсцисс, то и метод расходится;

- если, т.е. касательная к графику почти параллельна оси ординат, то и продвижения к корню не будет.

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

(11)

1.3 Решение нелинейного уравнения методом простых итераций

Нелинейное уравнение:

1/?(x-4)?^2 -0.5=0. (12)



Для приблизительной оценки корней нелинейного уравнения (12) построим график функции с помощью программы Matlab (рисунок 1.3).

Рисунок 1.3. График нелинейной функции y = f(x)

Из рисунка 1.3 видно, что график функции пересекает ось абсцисс в двух точках х1, расположенной на отрезке [2; 3], и х2, расположенной на отрезке [5; 5.5].

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

x=-8:0.1:8;

y=(x-4).^(-2) - 0.5;

figure(1);

plot (x, y, 'r');

axis([0 6 -1 4]);

xlabel('x');

ylabel('y');

grid on;

title ('метод простых итераций для F(x)=1/(x-4)^2-0.5');

eps=0.001;

kmax=20;

x=input ('введите приближение');

y=(x-4).^(-2) - 0.5;

if (y<0)

Листинг 1.1, лист 1 - Код программы, реализующей решение нелинейного уравнения методом простых итераций

x1=x-y;

D=x-x1;

k=0;

while (abs(D)>eps)&&(k<kmax)

x=x1;

y=(x-4).^(-2) - 0.5;

if (y>0)

x1=x+y;

end

if (y<0)

x1=x-y;

end

D=x-x1;

k=k+1;

disp(x1);

end

disp ('Корень:');

disp(x1);

disp ('Номер итерации:');

disp(k);

disp('Точность');

disp(D);

end

if (y>0)

x1=x+y;

D=x-x1;

k=0;

while (abs(D)>eps)&&(k<kmax)

x=x1;

y=(x-4).^(-2) - 0.5;

x1=x+y;

D=x-x1;

k=k+1;

disp(x1);

end

disp ('Корень:');

disp(x1);

disp ('Номер итерации:');

disp(k);

disp('Точность');

disp(D);

end

Листинг 1.1, лист 2



Рисунок 1.4. Результат работы программы

Результат работы программы представлен на рисунке 1.4. В данном случае итерационный процесс сходится; первый корень уравнения x1= 2.5856 найден с заданной точностью за 7 итераций, погрешность решения составляет 4.1928e-04; второй корень уравнения х2=5.4145 найден с заданной точностью за 5 итераций, погрешность решения составляет 5.9527е-04.

1.4 Решение нелинейного уравнения методом Ньютона

Нелинейное уравнение:

?exp (x+1) - 2.5-x???=0.? (13)

Для приблизительной оценки корней нелинейного уравнения (13) построим график функции с помощью программы MatLab (рисунок 1.5).



Рисунок 1.5. График нелинейной функции y = f(x)

График функции пересекает ось абсцисс в точках, находящихся на отрезке [0; - 2.5].

С помощью программы Matlab напишем код для решения нелинейного уравнения методом Ньютона.

x=-3:0.1:3;

y=exp (x+1) - 2.5-x;

figure(1);

plot (x, y, 'b');

axis([-3 2.5 -2 20]);

xlabel('x');

ylabel('y');

grid on;

title ('метод Ньютона для F(x)=exp (x+1) - 2.5-x');

eps=0.001;

kmax=30;

x=-3;

[F, F1]=fun1 (x);

[x1, D]=Newton (x, F, F1);

k=0;

while (abs(D)>eps)&&(k<kmax)

x=x1;

[F, F1]=fun1 (x);

[x1, D]=Newton (x, F, F1);

k=k+1;

disp(x1);

end

disp ('Корень:');

disp(x);

disp ('Номер итерации:');

disp(k);

Листинг 1.2, лист 1 - Код программы, реализующей решение нелинейного уравнения методом Ньютона

disp('Точность');

disp(D);

function [F, F1] = fun1 (x)

F=exp (x+1) - 2.5-x;%расчет ф-ции

if nargout>1

F1=exp (x+1) - 1;%1 производная

end

function [x1, D] = Newton (x, F, F1)

D=-F/F1;

x1=x+D;

end

Листинг 1.2, лист 2



Рисунок 1.6. Результат работы программы

Результат работы программы представлен на рисунке 1.6. Один из корней уравнения x= -2.1992 найден с заданной точностью за 2 итерации погрешность решения составляет 9.0035е-04.

Были рассмотрены два метода решения нелинейных уравнений: Ньютона и простых итераций. Метод простых итераций является простейшим для нахождения корней нелинейных уравнений. Он уступает в скорости сходимости методу Ньютона, который в свою очередь является модифицированным методом простых итераций.?



2. Интерполяция функций

Пусть функция f(x) задана таблицей своих значений x_i ?, y?_i на интервале [a; b]:

? y?_i=f (x_i), i=0,1,…, n, a?x_i?b. (14)

Задача интерполяции - найти функцию f(x), принимающую в точках x_i те же значения y_i.

Условие интерполяции:

F (x_i)=y_i, (15)

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

Если F(x) ищется только на отрезке [a; b] - то это задача интерполяции, а если за пределами первоначального отрезка, то это задача экстраполяции.

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

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

- аппроксимация - определение в явном виде параметров функции, описывающей распределение точек.

Задача нахождения интерполяционной функции F(x) имеет много решений, так как через заданные точки x_i ?, y?_(i) можно провести бесконечно много кривых, каждая из которых будет графиком функции, для которой выполнены все условия интерполяции. Для практики важен случай аппроксимации функции многочленами:



F(x)=P_m ?(x?_i)=a_0+a_1*x+a_2*x^2+ …+a_m*x^m, i=0,1,…, m, (16)

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

При построении одного многочлена для всего рассматриваемого интервала [a; b], для нахождения коэффициентов многочлена необходимо использовать все уравнения системы (16). Данная система содержит n+1 уравнение, следовательно, с ее помощью можно определить n+1 коэффициент. Поэтому максимальная степень интерполяционного многочлена m=n, и многочлен принимает вид:

? P?_n ?(x?_i)=a_0+a_1*x+a_2*x^2+ …+a_n*x^n, i=0,1,…, n. (17)

2.1 Локальная и глобальная интерполяция

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

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

Если функция f(x) интерполируется на отрезке [a; b] с помощью единого многочлена P_m (x) для всего отрезка, то такую интерполяцию называют глобальной. В случае локальной интерполяции на каждом интервале [x_i, x_(i+1)] строится отдельный интерполяционный полином невысокой степени.



2.2 Кусочно-линейная интерполяция

Простейшим и часто используемым видом локальной интерполяции является линейная (или кусочно-линейная) интерполяция. Она заключается в том, что узловые точки соединяются отрезками прямых (Рис. 3.1), то есть через каждые две точки (x_i, y_i) и (x_(i+1), y_(y+1)) проводится полином первой степени:

F(x)=a_0+a_1*x, x_(i-1)?x?x_i. (18)

Коэффициенты a_0 и a_1 разные на каждом интервале [x_i, x_(i+1)], и находятся из выполнения условий интерполяции на концах отрезка:

{f_(i-1)=a_0+a_1*x_(i-1);@f_i=a_0+a_1*x_i.)+ (19)

Из системы уравнений (19) можно найти коэффициенты:

a_0=f (x_(i-1)) - a_1*x_(i-1), a_1=(f (x_i) - f (x_(i-1)))/(x_i-x_(i-1)). (20)

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

Рисунок 2.1. Кусочно-линейная интерполяция



2.3 Кусочно-квадратичная интерполяция

В случае квадратичной интерполяции, для каждых трех узловых точек (x_(i-1), y_(y-1)), (x_i, y_y), (x_(i+1), y_(y+1)) строится уравнение параболы:

F(x)=a_0+a_1*x+a_2*x^2, x_(i-1)?x?x_(i+1). (21)

Здесь коэффициенты a_0 ?, a?_(1) и a_2 и разные на каждом интервале [x_(i-1), x_(i+1)] и определяются решением системы уравнений для условия прохождения параболы через три точки:

{(f_(i-1)=a_0+a_1*x_(i-1)+a_2*x_(i-1)^2@f_i=a_0+a_1*x_i+a_2*x_i^2@f_(i+1)=a_0+a_1*x_(i+1)+a_2*x_(i+1)^2)+ (22)

Из системы уравнений (22) можно найти коэффициенты:

a_0=f (x_(i-1)) - a_1*x_(i-1) - a_2*x_(i-1)^2,

a_1=(f (x_(i+1)) - f (x_(i-1)))/((x_(i+1) - x_(i-1))*(x_(i+1) - x_i)) - a_2*(x_i+x_(i-1)),

a_2=(f (x_(i+1)) - f (x_(i-1)))/((x_(i+1) - x_(i-1))*(x_(i+1) - x_i)) - (f (x_i) - f (x_(i-1)))/((x_i-x_(i-1))*(x_(i+1-) x_i)).

Рисунок 2.2. Кусочно-квадратичная интерполяция



2.4 Кусочно-квадратичная интерполяция функции по таблице значений

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

Таблица 2.1. Экспериментальные данные для интерполяции

0	Экспериментальные точки	Точки вычисления f(x)
1	x	-8	-6	-3,5	-3	-2,5	0	2	2,5	4	6,5	-10	-9	-5	-1	1,5	3	5	9
	y	-1	3	6,5	4	2	4	4,5	1	-2	1
2	x	-9,5	-6,5	-4	-2,5	-0,5	1,5	3	4,5	9,5		-8	-5	-3	-1,5	0,5	4	8
	y	5,5	-1	-4,5	-2	-5	-2,5	0	1,5	-1,5

С помощью программы Matlab напишем программу для интерполяции первого набора значений.

X=[-8 -6 -3.5 -3 -2.5 0.1 2 2.5 4 6.5];

F=[-1 3 6.5 4 2 4 4.5 1 -2 1];

g=[-10 -9 -5 -1 1.5 3 5 9];

i=1;

x=zeros();

e=zeros();

G=zeros();

E=zeros();

plot (X, F, 'or');

hold on;

plot (X, F);

hold on;

grid on;

n=10;

r=1;

H=func;

while (i<=8)

if (g(i)>=X(1)&&g(i)<=X(10)&&H(i)~=0)

j=H(i);

Листинг 2.1, лист 1 - Код программы, реализующей кусочно-квадратичную интерполяцию для первого набора данных

G(r)=g(i);

x0=X(j);

x1=X (j+1);

x2=X (j+2);

f0=F(j);

f1=F (j+1);

f2=F (j+2);

A=[x0^2 x0 1 f0; (x1)^2 x1 1 f1; (x2)^2 x2 1 f2];

p=1;

n=3;

k=1;

eps=0.001;

while (k<=n)

if (abs(A (k, k))>eps)

j=n+1;

while (j>=k)

A (k, j)=A (k, j)/A (k, k);

j=j-1;

end

p=k+1;

while (p<=n)

j=n+1;

while (j>=k)

A (p, j)=A (p, j) - A (p, k)*A (k, j);

j=j-1;

end

p=p+1;

end

k=k+1;

end

end

C=[A (1,1) A (1,2) A (1,3); A (2,1) A (2,2) A (2,3); A (3,1) A (3,2) A (3,3)];

d=[A (1,4); A (2,4); A (3,4)];

x(n)=d(n);

x (n-1)=d (n-1) - C (n-1, n)*x(n);

x (n-2)=d (n-2) - (C (n-2, n-1)*x (n-1)+C (n-2, n)*x(n));

a=x (n-2);

b=x (n-1);

c=x(n);

E=a*g(i)^2+b*g(i)+c;

r=r+1;

plot (g(i), E, '*r');

hold on;

end

i=i+1;

end

l=1;

u=1;

o=2;

m=zeros();

m(1)=G(u);

Листинг 2.1, лист 2

e(u)=G(u);

while (e(u)<G(5))

u=u+1;

e(u)=G(1)+0.01*(u-1);

if (e(u)~=0)

m(o)=e(u);

o=o+1;

end

end

i=1;

W=zeros();

H=func01 (m);

while (i<=1000)

if (H(i)>0)

j=H(i);

x0=X(j);

x1=X (j+1);

x2=X (j+2);

f0=F(j);

f1=F (j+1);

f2=F (j+2);

A=[x0^2 x0 1 f0; (x1)^2 x1 1 f1; (x2)^2 x2 1 f2];

p=1;

n=3;

k=1;

eps=0.0001;

while (k<=n)

if (abs(A (k, k))>eps)

j=n+1;

while (j>=k)

A (k, j)=A (k, j)/A (k, k);

j=j-1;

end

p=k+1;

while (p<=n)

j=n+1;

while (j>=k)

A (p, j)=A (p, j) - A (p, k)*A (k, j);

j=j-1;

end

p=p+1;

end

k=k+1;

end

end

C=[A (1,1) A (1,2) A (1,3); A (2,1) A (2,2) A (2,3); A (3,1) A (3,2) A (3,3)];

d=[A (1,4); A (2,4); A (3,4)];

x(n)=d(n);

x (n-1)=d (n-1) - C (n-1, n)*x(n);

x (n-2)=d (n-2) - (C (n-2, n-1)*x (n-1)+C (n-2, n)*x(n));

a=x (n-2);

Листинг 2.1, лист 3

b=x (n-1);

c=x(n);

r=r+1;

E=a*m(i)^2+b*m(i)+c;

W(i)=E;

end

i=i+1;

end

plot (m, W);

hold on;

function H = func

x0=[-8 -6 -3.5 -3 -2.5 0.1 2 2.5 4 6.5];

x=[-10 -9 -5 -1 1.5 3 5 9];

n=10, k=1;

H=zeros (1,8);

while (k<=8)

i=1;

while (i<n)

if (x(k)>x0 (i)&& x(k)<x0 (i+1))

H(k)=i-1;

end

i=i+1;

end

k=k+1;

end

end

function H = func01 (m)

x0=[-8 -6 -3.5 -3 -2.5 0.1 2 2.5 4 6.5];

n=10;

k=1;

H=zeros();

while (k<=1000)

i=1;

while (i<n)

if (m(k)>=x0 (i)&& m(k)<=x0 (i+1))

H(k)=i-1;

end

i=i+1;

end

k=k+1;

end

end

Листинг 2.1, лист 4

Графики функции, построенные с помощью программы Matlab, представлены на рисунке 2.3.

Рисунок 2.3. График функции для первого набора данных

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



Таблица 2.2. Значения функции в расчетных точках

	Точки вычисления f(x)
x	-10	-9	-5	-1	1.5	3	5	9
y			4.6	0,6154	4,4471	-1,25	-2

С помощью программы Matlab напишем программу для интерполяции второго набора значений.

X=[-9.5 -6.5 -4 -2.5 -0.5 1.5 3 4.5 9];

F=[5.5 -1 -4.5 -2 -5 -2.5 0.1 1.5 -1.5];

g=[-8 -5 -3 -1.5 0.5 4 8];

i=1;

x=zeros();

e=zeros();

G=zeros();

W=zeros();

plot (X, F, 'og');

hold on;

Листинг 2.2, лист 1 - Код программы, реализующей кусочно-квадратичную интерполяцию для второго набора данных

plot (X, F, 'g');

hold on;

grid on;

H=func2;

r=1;

while (i<=7)

if (g(i)>=X(1)&&g(i)<=X(9)&&H(i)~=0)

j=H(i);

G(r)=g(i);

x0=X(j);

x1=X (j+1);

x2=X (j+2);

f0=F(j);

f1=F (j+1);

f2=F (j+2);

A=[x0^2 x0 1 f0; (x1)^2 x1 1 f1; (x2)^2 x2 1 f2];

p=1;

n=3;

k=1;

eps=0.0001;

while (k<=n)

if (abs(A (k, k))>eps)

j=n+1;

while (j>=k)

A (k, j)=A (k, j)/A (k, k);

j=j-1;

end

p=k+1;

while (p<=n)

j=n+1;

while (j>=k)

A (p, j)=A (p, j) - A (p, k)*A (k, j);

j=j-1;

end

p=p+1;

end

k=k+1;

end

end

C=[A (1,1) A (1,2) A (1,3); A (2,1) A (2,2) A (2,3); A (3,1) A (3,2) A (3,3)];

d=[A (1,4); A (2,4); A (3,4)];

x(n)=d(n);

x (n-1)=d (n-1) - C (n-1, n)*x(n);

x (n-2)=d (n-2) - (C (n-2, n-1)*x (n-1)+C (n-2, n)*x(n));

a=x (n-2);

b=x (n-1);

c=x(n);

E=a*g(i)^2+b*g(i)+c;

g(i);

r=r+1;

Листинг 2.2, лист 2

plot (g(i), E, '*m');

hold on;

end

i=i+1;

end

l=1;

m=zeros();

u=1;

o=1;

e(u)=G(u);

while (e(u)<G(7))

u=u+1;

e(u)=G(1)+0.05*(u-1);

if (e(u)~=0)

m(o)=e(u);

o=o+1;

end

end

i=1;

H=func02 (m);

while (i<=319)

if (H(i)>0)

j=H(i);

x0=X(j);

x1=X (j+1);

x2=X (j+2);

f0=F(j);

f1=F (j+1);

f2=F (j+2);

A=[x0^2 x0 1 f0; (x1)^2 x1 1 f1; (x2)^2 x2 1 f2];

p=1;

n=3;

k=1;

eps=0.0001;

while (k<=n)

if (abs(A (k, k))>eps)

j=n+1;

while (j>=k)

A (k, j)=A (k, j)/A (k, k);

j=j-1;

end

p=k+1;

while (p<=n)

j=n+1;

while (j>=k)

A (p, j)=A (p, j) - A (p, k)*A (k, j);

j=j-1;

end

p=p+1;

end

k=k+1;

end

Листинг 2.2, лист 3

end

C=[A (1,1) A (1,2) A (1,3); A (2,1) A (2,2) A (2,3); A (3,1) A (3,2) A (3,3)];

d=[A (1,4); A (2,4); A (3,4)];

x(n)=d(n);

x (n-1)=d (n-1) - C (n-1, n)*x(n);

x (n-2)=d (n-2) - (C (n-2, n-1)*x (n-1)+C (n-2, n)*x(n));

a=x (n-2);

b=x (n-1);

c=x(n);

r=r+1;

E=a*m(i)^2+b*m(i)+c;

W(i)=E;

end

i=i+1;

end

plot (m, W);

hold on;

function H = func02 (m)

x0=[-9.5 -6.5 -4 -2.5 -0.5 1.5 3 4.5 9];

n=9;

k=1;

H=zeros();

while (k<=319)

i=1;

while (i<=n)

if (m(k)>=x0 (i)&& m(k)<=x0 (i+1))

if (i~=1)

H(k)=i-1;

else

H(k)=i;

end

end

i=i+1;

end

k=k+1;

end

end

function H = func2

x0=[-9.5 -6.5 -4 -2.5 -0.5 1.5 3 4.5 9];

x=[-8 -5 -3 -1.5 0.5 4 8];

n=9;

k=1;

H=zeros (1,7);

while (k<=7)

i=1;

while (i<n)

if (x(k)>x0 (i)&& x(k)<x0 (i+1))

if (i~=1)



Листинг 2.2, лист 4

H(k)=i-1;

else

H(k)=i;

end

end

i=i+1;

end

k=k+1;

end

end

Листинг 2.2, лист 5

Графики функции, построенные с помощью программы Matlab, представлены на рисунке 2.4.

Рисунок 2.4. График функции для второго набора данных



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

Таблица 2.3. Значения функции в расчетных точках

	Точки вычисления f(x)
x	-8	-5	-3	-1,5	0.5	4	8
y	1,9364	-3,3091	-3,2167	-2,5952	-4,4375	1,1667	0,1

2.5 Кубический сплайн

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

Таблица 2.4. Экспериментальные данные для интерполяции

0	Экспериментальные точки	Точки вычисления f(x)
1	x	-8	-6	-3,5	-3	-2,5	0	2	2,5	4	6,5	-10	-9	-5	-1	1,5	3	5	9
	y	-1	3	6,5	4	2	4	4,5	1	-2	1
2	x	-9,5	-6,5	-4	-2,5	-0,5	1,5	3	4,5	9,5		-8	-5	-3	-1,5	0,5	4	8
	y	5,5	-1	-4,5	-2	-5	-2,5	0	1,5	-1,5

Код программы в Matlab для построения графика неизвестной функции f(x) и вычисления значений функции в промежуточных точках для первого набора экспериментальных данных представлен ниже:

function SP (x0, y0);

n = length(x0) -1;

y = zeros (1, n-1);

h= zeros (1, n);

D = zeros (1, n-1);

A = zeros (1, n-1);

B = zeros (1, n-1);

C = zeros (1, n-1);

i = 1;

while (i<=n)

h(i)=x0 (i+1) - x0 (i);

i=i+1;

end

i=1;

while (i<=n-1)

D(i)=6*((y0 (i+2) - y0 (i+1))/h (i+1) - (y0 (i+1) - y0 (i))/h(i));

i=i+1;

end

A(1)=0;

i=2;

while (i<=n-1)

A(i)=h(i);

i=i+1;

end

i=1;

while (i<=n-1)

B(i)=2*(h(i)+h (i+1));

i=i+1;

end

i=1;

while (i<=n-2)

C(i)=h (i+1);

i=i+1;

end

Q=zeros (1, n);

R=zeros (1, n);

i=1;

while (i<=n-1)

Листинг 2.3, лист 1 - Код программы, реализующей интерполяцию кубическим сплайном для первого набора значений

Q (i+1)= - (C(i)/(B(i)+A(i)*Q(i)));

R (i+1)=(D(i) - A(i)*R(i))/(B(i)+A(i)*Q(i));

i=i+1;

end

M=zeros (1, n-1);

M (n-1)=R(n);

i=n-2;

while (i>=1)

M(i)=Q (i+1)*M (i+1)+R (i+1);

i=i-1;

end

M =[0, M, 0];

k=1;

while (k<n)

i=1;

while (i<=n)

l=x0 (i):0.1:x0 (i+1);

m=(1/(6*h(i)))*((M(i)*(x0 (i+1) - l).^3)+M (i+1)*(l-x0 (i)).^3)+(1/h(i))*((y0 (i) - ((M(i)*h(i)^2)/6))*(x0 (i+1) - l)+(y0 (i+1) - ((M (i+1)*h(i)^2)/6))*(l - (x0 (i))));

plot (l, m, 'b');

hold on;

i=i+1;

end

k=k+1;

end

plot (x0, y0,'r');

plot (x0, y0,'or');

xlabel('x0');

ylabel('y0');

grid on;

Листинг 2.3, лист 2

Графики функции, построенные с помощью программы Matlab, представлены на рисунке 2.5.

Рисунок 2.5. График функции для первого набора значений

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

Таблица 2.5. Значения функции в расчетных точках

	Точки вычисления f(x)
x	-10	-9	-5	-1	1,5	3	5	9
y	-3,8941	-2,4471	6,3306	1,2921	6,3429	-1,2670	-1,1245	4,8450

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



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

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

Методы их решения подразделяются на два класса:

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

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

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

Решить дифференциальное уравнение

(23)

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



3.1 Метод Эйлера

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

Пусть дано обыкновенное дифференциальное уравнение (ОДУ) с начальными условиями (задача Коши):

(24)

и удовлетворяются условия существования и единственности решения.

Требуется найти решение задачи Коши (24) на отрезке. Найдем решение в виде таблицы. Для этого разобьем отрезок на n равных частей и построим последовательность

,

где - шаг интегрирования.

Проинтегрируем исходное уравнение на отрезке:

Полученное соотношение можно переписать как

(25)



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

.

Подставляя полученный результат в формулу (25), получим основную расчетную формулу метода Эйлера:

(26)

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

(27)

Далее, определяя значение аргумента x по формуле, используя найденное значение y1 и полагая в формуле (26), вычисляем следующее приближенное значение интегральной кривой как

(28)

Поступая аналогичным образом при, определяем все остальные значения yk, в том числе последнее значение

(29)



которое соответствует значению аргумента.

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

Метод Эйлера может быть применен к решению систем дифференциальных уравнений.

Пусть задана система двух уравнений первого порядка:

(30)

с начальными условиями

.

Необходимо найти решение этой задачи Коши. Проводя аналогичные рассуждения, получаем расчетные формулы вида:

(31)

где h - шаг интегрирования.

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



3.2 Метод Эйлера-Коши

Метод Эйлера-Коши также относится к методам второго порядка и тоже требует двукратного вычисления функции f (x, y):

(y_(i+1)) Ю=y_i+h*f (x_i, y_i), (32)

y_(i+1)=y_i+h/2*(f (x_i, y_i)+f (x_(i+1), (y_(i+1))) Ю), i = 0, 1,…, n - 1. (33)

Метод Эйлера-Коши имеет погрешность порядка O(h2).

3.3 Решение ОДУ методом Эйлера

Решите ОДУ первого порядка (задачу Коши) методом Эйлера:

а) запрограммируйте метод Эйлера в программу Matlab и получите таблицу значений неизвестной функции y(x);

б) аппроксимируйте полученные данные любым способом и постройте приближенный график функции y(x) в программе Matlab.

Таблица 3.1. Задача Коши для решения методом Эйлера

	ОДУ и начальные условия	Отрезок	Шаг
0		от 0 до 2	0,2

Код метода в Matlab представлен ниже:

t0=0:0.2:2;

n=length(t0);

y0=zeros();

y0 (1)=4;

h=0.2;

t=t0;

i=1;

for t=0.2:0.2:2

y0 (i+1)=y0 (i)+h*(1/(t+2)+2*log (2*t+4)+2*cos (2*t) - 2*sin (2*t) - 2*y0 (i));

i=i+1;

end

SP3 (t0, y0);

y0

t0

function SP3 (t0, y0)

n = length(t0) -1;

h=ones (1, n);

D = ones (1, n-1);

A = ones (1, n-1);

B = ones (1, n-1);

C = ones (1, n-1);

i = 1;

while (i<=n)

h(i)=t0 (i+1) - t0 (i);

i=i+1;

end

i=1;

while (i<=n-1)

D(i)=6*((y0 (i+2) - y0 (i+1))/h (i+1) - (y0 (i+1) - y0 (i))/h(i));

i=i+1;

end

A(1)=0;

i=2;

while (i<=n-1)

A(i)=h(i);

i=i+1;

end

i=1;

while (i<=n-1)

B(i)=2*(h(i)+h (i+1));

i=i+1;

end

i=1;

while (i<=n-2)

C(i)=h (i+1);

i=i+1;

end

Q=zeros (1, n);

R=zeros (1, n);

Листинг 3.1, лист 1 - Код программы, реализующей решение ОДУ первого порядка методом Эйлера

i=1;

while (i<=n-1)

Q (i+1)= - (C(i)/(B(i)+A(i)*Q(i)));

R (i+1)=(D(i) - A(i)*R(i))/(B(i)+A(i)*Q(i));

i=i+1;

end

M=zeros (1, n-1);

M (n-1)=R(n);

i=n-2;

while (i>=1)

M(i)=Q (i+1)*M (i+1)+R (i+1);

i=i-1;

end

M = [0, M, 0];

k=1;

while (k<n)

i=1;

while (i<=n)

l=t0 (i):0.02:t0 (i+1);

m=(1/(6*h(i)))*((M(i)*(t0 (i+1) - l).^3)+M (i+1)*(l-t0 (i)).^3)+(1/h(i))*((y0 (i) - ((M(i)*h(i)^2)/6))*(t0 (i+1) - l)+(y0 (i+1) - ((M (i+1)*h(i)^2)/6))*(l - (t0 (i))));

plot (l, m, 'b');

hold on;

i=i+1;

end

k=k+1;

end

hold on;

t0

y0

plot (t0, y0,'r');

plot (t0, y0,'o');

xlabel('t');

ylabel('y');

grid on;

end

Листинг 3.1, лист 2

График неизвестной функции представлен на рисунке 3.1



Рисунок 3.1. График неизвестной функции y(x)

Результат работы программы представлен на рисунке 3.2.

Рисунок 3.2. Найденные значения функции для заданных точек

3.4 Решение ОДУ методом Эйлера-Коши

Решите ОДУ первого порядка (задачу Коши) методом Эйлера-Коши:

а) запрограммируйте метод Эйлера-Коши в программу Matlab и получите таблицу значений неизвестной функции y(x);

б) аппроксимируйте полученные данные любым способом и постройте приближенный график функции y(x) в программе Matlab.



Таблица 3.2. Задача Коши для решения методом Эйлера-Коши

	ОДУ и начальные условия	Отрезок	Шаг
0		от 0 до 4	0,4

Код метода в Matlab представлен ниже:

t0=0:0.4:4;

n=length(t0);

Листинг 3.2, лист 1 - Код программы, реализующей решение ОДУ методом Эйлера-Коши

y0=zeros (1, n);

y1=zeros (1, n);

y0 (1)=1;

h=0.4;

i=1;

t=t0;

for t=0:0.4:3.6

y1 (i+1)=y0 (i)+h*(y0 (i) - 2*sin(t) - cos (2*t)+2*cos(t)*(1-2*sin(t)));

t1=t+h;

y0 (i+1)=y0 (i)+(h/2)*(y0 (i) - 2*sin(t) - cos (2*t)+2*cos(t)*(1-2*sin(t))+(y1 (i+1) - 2*sin(t1) - cos (2*t1)+2*cos(t1)*(1-2*sin(t1))));

i=i+1;

end

t0

y0

lagr (t0, y0);

function lagr (t0, y0)

t=0:0.04:4;

n = length(t0);

n1 = length(t);

L=ones (1, n1);

k=1;

l=length(t);

while (k<=l)

S=0;

i=1;

while (i<=n)

P=1;

j=1;

while (j<=n)

C=t0 (i) - t0 (j);

if (i==j)

j=j+1;

else

P=P*(t(k) - t0 (j))/C;

j = j+1;

end

end

S=S+P*y0 (i);

i=i+1;

end

L(k)=S;

k=k+1;

end

hold on;

plot (t, L, 'b');

xlabel('t');

ylabel('y');

plot (t0, y0,'r');

end

Листинг 3.2, лист 2

График неизвестной функции представлен на рисунке 3.3.

Рисунок 3.3. График неизвестной функции y(x)

Результат работы программы представлен на рисунке 3.4.

Рисунок 3.4. Результат работы программы

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

Метод Эйлера-Коши имеет те же недостатки, что и метод Эйлера, однако он обладает более высокой точностью, так как его погрешность пропорциональна шагу h во второй степени.



Заключение

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

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

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

Подчеркнем важные отличия численных методов от аналитических.

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

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

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



Библиографический список

1. Джон Г. Мэтьюз, Куртис Д. Финк. Численные методы. Использование Matlab. M.: Вильямс, 2001. - 716 с.

2. Вержбицкий В.М. Численные методы. Математический анализ и обыкновенные дифференциальные уравнения. М.: Оникс 21 век, 2005. - 400 с.

3. Интерполяция [Электронный ресурс] Режим доступа: http://aco.ifmo.ru/el_books/numerical_methods/lectures/glava3.html

4. Метод Эйлера-Коши [Электронный ресурс] Режим доступа: http://vm.psati.ru/online-vmath/index.php? page=27

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