Численные методы в среде MATLAB. Интерполяция функций одной и нескольких переменных.

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

На практике часто приходится использовать табличные данные вида {xi, yi} для приближенного вычисления у при любом значении аргумента х. Этой цели служит задача о приближении (аппроксимации) функций: данную функцию f(x) требуется приближенно заменить некоторой функцией g(x) так, чтобы отклонение g(x) от f(x) в заданной области было наименьшим. Если приближение строится на заданном дискретном множестве точек {хi}, то аппроксимация называется точечной. К ней относятся интерполирование, среднеквадратичное приближение и др. Для решения подобной задачи довольно часто используют сплайн-интерполяцию. Один из наиболее распространенных вариантов – интерполяция кубическими сплайнами. Кроме того существуют квадратичные и линейные сплайны. В MatLab интерполяция табличной функции одной переменной может быть произведена с помощью функции interp1(x, y, xi, method), где x, y – массивы одинакового размера 1хn, задающие исходную таблицу значений аргумента и функции, xi – массив 1xm, содержащий координаты точек, в которых необходимо найти значения функции, method – строковый параметр, определяющий способ интерполяции: 'nearest' кусочно-постоянная интерполяция 'linear' - кусочно-линейная интерполяция 'spline' - кусочно непрерывная интерполяция кубическим сплайном (по умолчанию) 'cubic' - кусочно-непрерывная интерполяция кубическим полиномом

Например:

> > x = 0: 10; y = sin(x); xi = 0:.25: 10;

> > yi = interp1(x, y, xi);

> > plot(x, y, 'o', xi, yi)

Коэффициенты сплайна в MatLab можно получить с помощью функции p = spline(x, y). При этом, p –является структурой с полями: form - строковый параметр (тип структуры), breaks - координаты узлов сплайна, coefs – искомые коэффициенты сплайна (матрица), pieces - количество интервалов, order – порядок, dim - размерность.

В Scilab для построения линейной интерполяции служит функция interpln(z, x), где z - матрица исходных данных, x – массив с координатами интересующих нас точек. Функция возвращает значения линейного сплайна в точках x. Например:

> > x=[132 140 150 162 170]; y=[330 350 385 425 450];

> > z=[x; y]; t=132: 5: 252;

> > ptd=interpln(z, t);

> > plot(t, ptd);

Построение кубического сплайна в SciLab состоит из двух этапов: вначале вычисляются коэффициенты сплайна с помощью функции d=splin(x, y), и только затем рассчитываются значения интерполяционного полинома в точках с абсциссами t с помощью функции у=interp(t, x, y, d). Эти функции имеют следующие параметры: x – строго возрастающий вектор, состоящий минимум из двухкомпонент, y – вектор того же размера, что и x, d – коэффициенты кубического сплайна, t – вектор абсцисс, y – вектор ординат, являющихся значениями кубического сплайна в точках t.

 

Численные методы в среде MATLAB. Решение алгебраических уравнений и систем уравнений.

Решение алгебраических уравнений и систем– solve

Для решения систем алгебраических уравнений и одиночных уравнений служит функция solve:

? solve(expr1, expr2,..., exprN, var1, var2,..., varN) – возвращает значения переменных varI, при которых соблюдаются равенства, заданные выражениями exprI. Если в выражениях не используются знаки равенства, то полагается exprI=0.

Результат решения возможен в следующих формах:

? для одного уравнения и одной переменной решение возвращается в виде одномерного или многомерного массива ячеек;

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

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

Решить уравнение

x3-1=0.

Решение:

> > syms x

> > y=x^3-1;

> > S=solve(y, x)

S =

[ 1]

[ -1/2+1/2*i*3^(1/2)]

[ -1/2-1/2*i*3^(1/2)]

Получили три корня x1=1, x2=, x3=, которые хранятся соответственно в элементах S(1), S(2), S(3) массива S. Подставим теперь их в этом порядке вместо x в выражение y=x3-1.

> > subs(y, x, S)

ans =

[ 0]

[ (-1/2+1/2*i*3^(1/2))^3-1]

[ (-1/2-1/2*i*3^(1/2))^3-1]

> > [m, n]=simple(ans)

m =

[ 0]

[ 0]

[ 0]

Выражение y=x3-1 принимает значение 0 при подстановке любого из корней, т.е. x1=1, x2=, x3=, являются точными корнями уравнения x3-1=0.

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

x5+x3+1=0.

Решение:

> > syms x

> > y=x^5+x^3+1;

> > S=solve(y, x);

> > vpa(S, 5)

ans =

[ -.83762]

[ -.21785-1.1670*i]

[ -.21785+1.1670*i]

[.63666-.66470*i]

[.63666+.66470*i]

Подставим теперь вектор S из 5-ти полученных значений корней вместо x в выражение

y=x5+x3+1.

> > subs(y, x, S)

ans =

[ 0]

[.9e-31+.1e-30*i]

[.9e-31-.1e-30*i]

[.1e-31+.1e-31*i]

[.1e-31-.1e-31*i]

Как видим, практически y=0, т.е. приближенные значения корней уравнения x5+x3+1=0 найдены с высокой степенью точности. Корни многочлена y=x5+x3+1 можно найти с помощью функции roots:, рассмотренной в

> > S=roots([1 0 1 0 0 1])

S =

-0.21785 + 1.167i

-0.21785 - 1.167i

0.63666 + 0.6647i

0.63666 - 0.6647i

-0.83762

Функция solve позволяет решать уравнения, представленные в аналитическом виде.

Решить уравнение

ax2+bx+c=0.

Решение:

> > S=solve('a*x^2+b*x+c=0', x)

S =

[ 1/2/a*(-b+(b^2-4*a*c)^(1/2))]

[ 1/2/a*(-b-(b^2-4*a*c)^(1/2))]

Это известные известные выражения корней x1, 2 = квадратного уравнения ax2+bx+c=0 через его коэффициенты a, b, c (обратная теорема Виета). Точно также можно выразить с помощью радикалов решения кубического уравнения ax3+bx2+cx+d=0, хотя эти выражения достаточно сложные.

Решение любого трансцендентного уравнения, в том числе и тригонометрического, достаточно сложная и серьезная проблема для систем аналитических вычислений. Иногда простое трансцендентное уравнение может и не решаться в MATLAB. Когда MATLAB не может найти ни одного решения, то функция solve возвращает пустую последовательность[ empty sym ]. Это означает, что либо решения не существует, либо MATLAB не удалось его найти. В таких случаях, возможно, системе следует помочь, сделав какие то нестандартные преобразования уравнения или системы, приведя их к другому виду.

Функция fzero позволяет найти все вещественные корни уравнения f(x)=0. Решить трансцендентное уравнение

ln(4-2x)+x2-2=0.

Решение:

> > syms x

> > Y=log(4-2*x)+x^2-2;

> > S=solve(Y)

S =

1.27740129788842

> > subs(Y, S)

ans =

-.1e-14

Итак, найден корень x1=1, 27740129788842. Преобразуем уравнение к виду

e2? x? +2x-4=0.

Решение:

> > syms x

> > Y=exp(2-x^2)+2*x-4

Y =

exp(-x^2+2)-4+2*x

> > S=solve(Y)

S =

1.90009059438035

> > subs(Y, S)

ans =

.7e-14

Найден второй корень x2=1, 90009059438035.

Функция fzero позволяет найти все вещественные корни уравнения f(x)=0. Для этого небходимо знать области изоляции корней. Для их определения достаточно с помощью функции ezplot (реализации этой функции рассмотрены далее) построить график функции y=f(x) (рис. ) и визуально определить области изоляции всех корней.

> > syms x

> > Y=log(4-2*x)+x^2-2;

> > ezplot(Y, [-2, 3]), grid

Из рис. видно, что уравнение

ln(4-2x)+x2-2=0

имеет три вещественных корня, а областями изоляции корней могут быть: [1; 1, 5], [1, 5; 1, 99], [-1; -0, 5]. Первых два корня определила функция solve. Найдем оставшийся корень помощью функции fzero.

> > fzero('log(4-2*x)+x^2-2', [-1, -0.5], 0.1e-15)

ans =

-0.59451826012754

disp(subs(Y, ans))

Найден третий корень x3=-0, 59451826012754.

Существенный недостаток функции fzero состоит в том, что она не определяет комплексных корней уравнения. Совместно функции solve и fzero могут найти все вещественные и комплексные корни.

Прейдем теперь к системам уравнений.

Решить систему уравнений

> > syms x y

> > Y1=x+y-3;

> > Y2=x*y^2-4;

> > S=solve(Y1, Y2, x, y)

S =

x: [3x1 sym]

y: [3x1 sym]

disp([S.x S.y])

[ 4, -1]

[ 1, 2]

[ 1, 2]

Получили два решения (x1 y1)=(4; -1) и (x2, y2)=(1; 2) (второе кратное 2), причем (x1; y1) хранится в [S.x(1) S.y(1)], а (x2; y2) – в [S.x(2) S.y(2)]:

> > disp([S.x(1) S.y(1)])

[ 4, -1]

> > disp([S.x(2) S.y(2)])

[ 1, 2]

Подставим в выражения Y1=x+y-3 и Y2=xy2-4 вначале первое решение, а затем второе.

> > disp(subs([Y1 Y2], [x y], [S.x(1) S.y(1)])

[ 0, 0]

> > disp(subs([Y1 Y2], [x y], [S.x(2) S.y(2)])

[ 0, 0]

Видим, что решения точные, т.к. выражения Y1 и Y2 при их подстановке обратились в 0.

Функция solve позволяет решать системы уравнений, представленные в аналитическом виде.

Решить систему уравнений

Решение:

> > syms a b c x y z

> > Y1=(a+b)/(x+y)+(b+c)/(y+z)-(c+a)/(z+x)-1;

> > Y2=(a+b)/(x+y)-(b+c)/(y+z)+(c+a)/(z+x)-1;

> > Y3=-(a+b)/(x+y)+(b+c)/(y+z)+(c+a)/(z+x)-1;

> > S=solve(Y1, Y2, Y3, x, y, z)

S =

x: [1x1 sym]

y: [1x1 sym]

z: [1x1 sym]

> > disp([S.x S.y S.z])

[ a, b, c]

> > subs([Y1 Y2 Y3], [x y z], [S.x S.y S.z])

ans =

[ 0, 0, (-a-b)/(a+b)+1]

> > disp(simplify(ans))

[ 0, 0, 0]

Итак, решение системы (x; y; z)=(a; b; c), в чем мы убедились подстановкой.

Приведем пример решения системы трансцендентных уравнений

Решение:

> > syms x y

> > Y1=3^y*9^x-81;

> > Y2=log10((y+x)^2)-log10(x)-2*log10(3);

> > S=solve(Y1, Y2, x, y)

S =

x: [4x1 sym]

y: [4x1 sym]S =

> > R=[S.x S.y];

> > disp(vpa(R, 10))

[ 16.00000002, -28.00000004]

[ 16.00000002, -3.999999992]

[ 1.000000000, -3.999999996]

[ 1.000000000, 1.999999996]

 

⇐ Предыдущая1234567Следующая ⇒

Поделиться: