Полиномиальная интерполяция

Материал из Xgu.ru

Перейти к: навигация, поиск
stub.png
Данная страница находится в разработке.
Эта страница ещё не закончена. Информация, представленная здесь, может оказаться неполной или неверной.

Если вы считаете, что её стоило бы доработать как можно быстрее, пожалуйста, скажите об этом.

Автор: Владимир Кореньков
Правильная ссылка: http://xgu.ru/wiki/octave/interpolation/polinom

Содержание


Icon-caution.gif

Внимание! Ниже представлен возможно не самый лучший способ решения задачи. Если у Вас есть замечания -- обязательно сообщите они них на данной страничке либо автору.

[править] Постановка задачи

Имеется уравнение вида:

y=\frac{82000}{40x^{0.5}}+\frac{923}{40 x^{0.073}+1}

Определить: обратную зависимость, т.е., функцию вида x = f(y)

Так как показатели степеней являются дробными, то очевидно, что системы символьных вычислений (типа Maxima, Axiom и пр.) с подобной задачей не справятся. Следовательно, одним из наиболее простых решений будет следующее:

  1. Табулировать функцию, т.е. задавшись некоторым диапазоном [x1,x2] вычислить значения [y1,y2].
  2. Интерполировать данные [xi = f(yi)] другой функцией.

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

  • Интерполировать табличные данные сплайнами. Этот способ хорош тем, что позволяет описать практически любые кривые. Однако полином, полученный в результате вычисления, займет при написании несколько страниц (если количество начальных значений, как в рассматриваемом случае, >1000).
  • Если заранее известна зависимость, то используя метод наименьших квадратов, можно вычислить коэффициенты данной зависимости. Например, судя по графику функции (см. рисунок ниже), можно сделать вывод о том, что кривая описывается уравнением вида y=a \cdot e^{bx}, в котором следует подобрать неизвестные a и b.
  • Использовать полиномиальную интерполяцию

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

[править] Решение

Пример 1. Использование встроенной функции "polyfit"

# Файл SolveEquations.m

# Определяем диапазон начальных значений: от ста до миллиона с шагом 10000
xx = 10^3:10^4:10^6;
# Вычисляем значение функции
yy = 82000./((40.*xx).^0.5).+923./((40.*xx).^0.073.+1) ;

# Меняем столбцы местами
x = yy;
y = xx;

# Методом подбора (повторного вызова данного скрипта)
# определяем степень полинома равной 8
n = 8;

# Вычисляем коэффициенты "а" полинома
# y = a(1)*x^8 + a(2)*x^7 + a(3)*x^6 + ... + a(8)*x + a(9)
a = polyfit(x,y,n);


# Переопределяем функцию y = Sum(a_{i} * x^{i+n+1}), для i=1..n+1
function result=fun(x,a)
result=0;
n = length(a)-1;
for (i=1:n+1)
result = result .+ a(i) .* x .^(n .+ 1 - i);
endfor
endfunction


# Вычисляем точки интерполяционного полинома для
# графического сравнения с исходными данными
z = fun(x,a);

# Строим графики
plot(x,y,"r",x,z,"b*");hold on; grid on;

# Формируем строку -- уравнение полинома
i = 0:n;
s = "";
for (i=1:n+1)
if (i==n+1)
s = strcat(s, num2str(a(i)) );
else
s = strcat(s, num2str(a(i)), "*x^", num2str(n .+ 1 - i), "+" );
endif
endfor

# Выводим уравнение на печать
s

Пример использования

> addpath ...путь к файлу
> SolveEquations
> s = 4.6859e-11*x^8+-1.3027e-07*x^7+0.00015482*x^6+-0.10323*x^5+42.389
  • x^4+-11005*x^3+1.7678e+06*x^2+-1.6096e+08*x^1+6.3706e+09
Interpolation.png
Octave GNU Octave

Инсталляция | Синтаксис языка | Командная строка
Скрипты | Функции | Регулярные выражения | Массивы | Графики | Ввод/вывод данных
Распределенные вычисления | Численные методы | Сплайны