Полиномиальная интерполяция
Материал из Xgu.ru
Автор: Владимир Кореньков
Правильная ссылка: http://xgu.ru/wiki/octave/interpolation/polinom
Содержание |
|
|
Внимание! Ниже представлен возможно не самый лучший способ решения задачи. Если у Вас есть замечания -- обязательно сообщите они них на данной страничке либо автору. |
[править] Постановка задачи
Имеется уравнение вида:
Определить: обратную зависимость, т.е., функцию вида x = f(y)
Так как показатели степеней являются дробными, то очевидно, что системы символьных вычислений (типа Maxima, Axiom и пр.) с подобной задачей не справятся. Следовательно, одним из наиболее простых решений будет следующее:
- Табулировать функцию, т.е. задавшись некоторым диапазоном [x1,x2] вычислить значения [y1,y2].
- Интерполировать данные [xi = f(yi)] другой функцией.
Если первый этап не вызывает никаких вопросов, то второй этап можно выполнять несколькими путями:
- Интерполировать табличные данные сплайнами. Этот способ хорош тем, что позволяет описать практически любые кривые. Однако полином, полученный в результате вычисления, займет при написании несколько страниц (если количество начальных значений, как в рассматриваемом случае, >1000).
- Если заранее известна зависимость, то используя метод наименьших квадратов, можно вычислить коэффициенты данной зависимости. Например, судя по графику функции (см. рисунок ниже), можно сделать вывод о том, что кривая описывается уравнением вида
, в котором следует подобрать неизвестные 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)
endfunction
- result=0;
- n = length(a)-1;
- for (i=1:n+1)
- result = result .+ a(i) .* x .^(n .+ 1 - i);
- endfor
# Вычисляем точки интерполяционного полинома для
# графического сравнения с исходными данными
z = fun(x,a);
# Строим графики
plot(x,y,"r",x,z,"b*");hold on; grid on;
# Формируем строку -- уравнение полинома
i = 0:n;
s = "";
for (i=1:n+1)
endfor
- if (i==n+1)
- s = strcat(s, num2str(a(i)) );
- else
- s = strcat(s, num2str(a(i)), "*x^", num2str(n .+ 1 - i), "+" );
- endif
# Выводим уравнение на печать
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
| GNU Octave Инсталляция |
Синтаксис языка |
Командная строка
|
|---|
