Вычислительные методы в инженерных расчетах

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

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

Решение, полученное численным методом, обычно является приближенным, т. е. содержит некоторую погрешность. Источниками погрешности приближенного решения являются: 1) несоответствие математической задачи изучаемому реальному явлению; 2) погрешность исходных данных; 3) погрешность метода решения; 4) погрешности округлений в арифметических и других действиях над числами.

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

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

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

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

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

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

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

Оделить корни уравнения графически.

Уточнить один из них (точность 10) несколькими методами:

а) методом половинного деления б) методом хорд в) методом Ньютона г) используя стандартную формулу Matlab.

1.

2.

3.

1.Численные решения нелинейных уравнений Определена f (x) на множестве gR. Найти f (x)=0 (1).

Методы решения (1) бывают прямыми и итерационными.

Задача поиска разбивается на два этапа:

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

Задания начального значения корня и его дальнейшее уточнение до достижения заданной точности.

Получение очередного значения X называется к-ой итерацией.

Если {X}C, то говорят что итерационный процесс сходиться.

Основные методы уточнения корней:

Метод половинного деления.

Метод хорд (секущих).

Метод Ньютона (касательных).

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

1.1 Отделения корней уравнения графически

1. Программа

x=0.2:0.01:1;

y=2.71.^(-2*x)-2*x+1;

plot (x, y);

grid on

pause

x1=1.9:0.1:2.1;

y1= sin (x1+(pi/3))-x½;

plot (x1,y1);

grid on

pause

x2=0.5:0.01:1.5;

y2=3*x2.^4+8*x2.^3+6*x2.^2+10;

plot (x2,y2);

grid on

2. Графики Рис. 1 — График зависимости y (x)

Рис. 2 — График зависимости y1(x1)

Рис. 3 — График зависимости y2(x2)

1.2 Уточнение корней (точность 10)

1.2.1. Метод половинного деления

1) f (x) — непрерывна и пройден этап локализации.

2) f (a)f (b)<0

Рис. 4

Если f (a)*f (x)<0, то [a, b]=[a, x];

Если f (b)*f (x)<0, то [a, b]=[x, b].

На n-ом шаге отрезок уменьшается в 2 раз и продолжаем пока:

.

1. Программа

function y=mode (x)

y= 2.71.^(-2*x)-2*x+1;

a=0.5;

b=1;

eps=10^(-6);

L=b-a;

m=1;

while L>eps;

x=(a+b)/2;

c (m)=x;

if mode (x)*mode (a)<0;

b=x;

else

a=x;

end

m=m+1;

L=b-a;

end

disp ('число шагов дл уточнени корн')

disp (m)

disp ('последовательность приблежений')

disp (c)

disp ('корень уровнени')

disp (x)

plot (c)

Рис. 5 — Решение Рис. 6 — График последовательности приближений с (m)

2. Программа

function y1=mode2(x1)

y1=sin (x1+(pi/3))-x½;

a=1.9;

b=2.1;

eps=10^(-6);

L=b-a;

m=1;

while L>eps;

x1=(a+b)/2;

c (m)=x1;

if mode2(x1)*mode2(a)<0;

b=x1;

else

a=x1;

end

m=m+1;

L=b-a;

end

disp ('число шагов дл уточнени корн')

disp (m)

disp ('последовательность приблежений')

disp (c)

disp ('корень уровнени')

disp (x1)

plot (c)

Рис. 7 — Решение Рис. 8 — График последовательности приближений c (m)

3. Программа

function y2=mode3(x2)

y2=3*x2.^4+8*x2.^3+6*x2.^2+10;

a=0.5;

b=1.5;

eps=10^(-6);

L=b-a;

m=1;

while L>eps;

x2=(a+b)/2;

c (m)=x2;

if mode3(x2)*mode3(a)<0;

b=x2;

else

a=x2;

end

m=m+1;

L=b-a;

end

disp ('число шагов дл уточнени корн')

disp (m)

disp ('последовательность приблежений')

disp (c)

disp ('корень уровнени')

disp (x2)

plot (c)

Рис. 9 — Решение Рис. 10 — График последовательности приближений c (m)

1.3 Метод хорд

1) f (x) — непрерывна на [a, b]

2) f (a)*f (b)<0

f (a)>0, f (b)<0

Рис. 11

Пологая у=0

)

Если f (a)*f (x)<0, то [a, b]=[a, x];

Если f (b)*f (x)<0, то [a, b]=[x, b].

Алгоритм в ряде случаев сходиться.

1. Программа

clc

a=0.5;

b=2;

eps=10^(-6);

x=(a-(b-a)/(mode (b)-mode (a))*mode (a));

m=1;

while abs (mode (x))>eps;

c (m)=x;

if mode (x)*mode (a)<0;

b=x;

else

a=x;

end

m=m+1;

x=(a-(b-a)/(mode (b)-mode (a))*mode (a));

end

disp ('число шагов дл уточнени корн')

disp (m)

disp ('последовательность приблежений')

disp (c)

disp ('корень уровнени')

disp (x)

plot (c)

Рис. 12 — Решение Рис. 13 — График последовательности приближений c (m)

2. Программа

a=1.9;

b=2.1;

eps=10^(-6);

L=b-a;

m=1;

while L>eps;

x1=(a+b)/2;

c (m)=x1;

if mode2(x1)*mode2(a)<0;

b=x1;

else

a=x1;

end

m=m+1;

L=b-a;

end

disp ('число шагов дл уточнени корн')

disp (m)

disp ('последовательность приблежений')

disp (c)

disp ('корень уровнени')

disp (x1)

plot (c)

Рис. 14 — Решение Рис. 15 — График последовательности приближений c (m)

3. Программа

clc

a=0.5;

b=1.5;

eps=10^(-6);

L=b-a;

m=1;

while L>eps;

x2=(a+b)/2;

c (m)=x2;

if mode3(x2)*mode3(a)<0;

b=x2;

else

a=x2;

end

m=m+1;

L=b-a;

end

disp ('число шагов дл уточнени корн')

disp (m)

disp ('последовательность приблежений')

disp (c)

disp ('корень уровнения')

disp (x2)

plot (c)

Рис. 16 — Решение Рис. 17 — График последовательности приближений c (m)

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

Функция f (x) дифференцирована в с

f'(x) не должна менять свой знак на промежутке от [c, x]

Метод Ньютона также называют методом касательных. Рассмотрим в точке х касательную к кривой y=f (x), задаваемую уравнением

.

Положив Y=0, находим точку x пересечение касательной с осью абсцисс:

Построив касательную в точке x, получаем Рис. 18

по аналогичной формуле точку x пересечения этой касательной с осью x и т. д.:

1. Программа

a=0.5;

b=1;

eps=10^(-6);

x=(a-(b-a)/(mode (b)-mode (a))*mode (a));

m=1;

while abs (mode (x))>eps;

x=x-mode (x)/modew (x);

c (m)=x;

m=m+1;

end

disp ('число шагов дл уточнени корн')

disp (m)

disp ('последовательность приблежений')

disp (c)

disp ('корень уровнени')

disp (x)

plot (c)

Рис. 19 — Решение Рис. 20 — График последовательности приближений с (m)

Программа

function v1=modew2(x1)

v1=cos (x1+(pi/3))+½;

a=1.9;

b=2.1;

eps=10^(-6);

x=(a-(b-a)/(mode2(b)-mode2(a))*mode2(a));

m=1;

while abs (mode2(x1))>eps;

x=x-mode2(x1)/modew2(x1);

c (m)=x;

m=m+1;

end

disp ('число шагов дл уточнени корн')

disp (m)

disp ('последовательность приблежений')

disp (c)

disp ('корень уровнени')

disp (x)

plot (c)

Рис. 21 — Решение Рис. 22 — График последовательности приближений c (m)

Программа

function v2=modew3(x2)

v2=2*(x2−1)+2.^x2*log (2);

a=0.5;

b=1.5;

eps=10^(-6);

x=(a-(b-a)/(mode3(b)-mode3(a))*mode3(a));

m=1;

while abs (mode3(x))>eps;

x=x-mode3(x)/modew3(x);

c (m)=x;

m=m+1;

end

disp ('число шагов дл уточнени корн')

disp (m)

disp ('последовательность приблежений')

disp (c)

disp ('корень уровнени')

disp (x)

plot (c)

Рис. 23 — Решение Рис. 24 — График последовательности приближений c (m)

2. Решение СЛАУ

Написать систему линейных уравнений 5Ч5

Проверить ее обусловленность Вычислить определитель матрицы Найти обратную матрицу Вычислить число обусловленностей в 3-х нормах Решить систему методом Зейделя.

2.1 Метод решения СЛАУ

Метод решения СЛАУ можно разделить на 2 класса:

1) прямые (точные)

2) итерационные.

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

1. Метод Крамера

2. Метод Гаусса

3. Метод прогонки

4. Метод простых итераций (метод Зейделя) Число обусловленности матрицы.

— число обусловленности матрицы А.

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

А*А=Е

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

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

Замечание

Aсимметрична

При решение СЛАУ полезно вычислить предварительно число обусловленностей и только потом подходить к решению.

2.2 Нормы векторов и матриц

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

Нормой Х=(х, х,…, х)' в nмерном вещественном пространстве R называется не отрицательная величина вычисляемая как функция компоненты вектора и обладающая следующими свойствами:

а) ()

б), для С) в)

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

а)

б) в)

г)

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

— симметрическая матрица

— вещественные

— максимальное значение модуля собственных значений матрицы Норма (3)-спектральная норма:

— если матрица A сама симметрическая () то, тогда =

— если имеем собственные значения

AX=X

Норма (3) для матриц оказывается равной:

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

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

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

Всегда можно обеспечить его сходимость. Также как и в методе простых итераций строиться эквивалентная система:

Xначальное значение;

— вектор решения на k-ой итерации.

B — нижняя треугольная матрица полученная из матрицы с диагональными элементами, равные 0.

С — верхний треугольная матрица полученная из матрицы

СЛАУ называется нормальной если ее матрица системы симметрическая т. е. A=A

(Ах, х)>0

линейный уравнение корень матрица Критерий Сильвестра:

Если, А — нормальная, то а>0 i=1,n.

Из любой не вырожденной матрицы можно сделать нормальную просто путем транспонирования Т-ма: метод Зейделя всегда сходится для нормальных СЛАУ к единственному решению при любом выборе начального приближения.

с>0 это позволяет привести систему уравнения к виду решения итерационного процесса СЛАУ Программа:

function [z1,z2]=zeidel2(A, B, eps)

N=size (A, 1);

C=A'*A;

D=A'*B;

for i=1:N

D1(i)=D (i)/C (i, i);

end

%D1=D1'

d1=D1;

for i=1:N

for j=1:N

if i==j

C1(i, j)=0;

else

C1(i, j)=-C (i, j)/C (i, i);

end

end

end

R1=d1;

K=0;

while K==0

for i=1:N

v=C1(i, 1: N);

a=dot (v', d1);

d1(i)=a+D1(i);

R2=d1;

S=max (abs (R2-R1));

if S

z1=d1';

z2=S;

K=K+1;

end;

R1=R2;

end

end

A=[1 -2 3 -4 5

7 8 -9 1 -2

4 -5 6 -7 8

1 2 -3 4 -5

7 -8 9 -1 2]

B=[6; 3; 9; 6; 3]

f=det (A)

A1=inv (A)

x=AB

slave=abs (A)

n1=max (sum (slave'))

slave1=abs (A1)

n11=max (sum (slave1'))

m1=n1*n11

n2=max (sum (slave))

n22=max (sum (slave1))

m2=n2*n22

r1=A*A'

Lambda1=eig (r1)

a1=max (abs (Lambda1))

n3=sqrt (a1)

r2=A1*A1'

Lambda2=eig (r2)

a2=max (abs (Lambda2))

n33=sqrt (a2)

m3=n3*n33

eps=10^(-6);

[z1,z2]=zeidel2(A, B, eps)

Решение:

A =

1 -2 3 -4 5

7 8 -9 1 -2

4 -5 6 -7 8

1 2 -3 4 -5

7 -8 9 -1 2

B =

m1 = 1.5441e+017

m2 = 2.7022e+017

m3 = 1.4352e+017

z1 = 0.5403

— 3.5372

— 2.7791

— 6.3206

— 4.8029

z2 = 9.9994e-007