Нахожу многочлен Лагранжа. На выходе получаю выражение вида:
L(x) = - a0(x - x1)(x - x2)(x - x3)(x - x4)(x - x5) + a1(x - x0)(x - x2)(x - x3)(x
- x4)(x - x5) - ... - a5(x - x0)(x - x1)(x - x2)(x - x3)(x - x4)
Необходимо привести к каноническому виду:
L(x) = A0 * x5 + ... + A4 * x + A5
Необходимо реализовать это программно.
В связи с чем вопрос: есть ли где готовые реализации подобного кода? На любом языке.
Мне бы посмотреть разобраться.
Уверен, что задача легко решаема, но буду рад любой реализации. Преподаватель внезапно
решил, что ручной вариант его не устраивает и надо программно, ещё и сроки небольшие,
а голова и без того забита вовсю курсовой.
UPD: Пролистайте до конца - там есть моя версия готовой программы.
Или нажмите сюда.
Ответы
Ответ 1
Для .Net существует библиотека под названием MathNet.Symbolics. Конечно, её возможности
намного меньше, чем у SymPy, но с вашей задачей тоже справится без проблем. MathNet.Symbolics
написана на F#, но вы также можете её использовать на любом другом языке из семейства .Net.
К слову, если вы предпочитаете платформу .NET, то вы можете работать с SymPy при
помощи IronPython
Мне не очень нравится идея вручную задавать выражение, поэтому для удобства я написал
простую функцию для его генерации
let test n =
let y = sprintf "A%i"
let l i n = seq {
for j in 0..n - 1 do
if j <> i then
yield sprintf "(x - x%i)" j
}
let p i =
let y' = y i
l i n
|> String.concat "*"
|> sprintf "%s*%s" y'
Array.init n p
|> String.concat "+"
open MathNet.Symbolics
let n = 4
let some = test n
printfn "%s" some
//A0*(x - x1)*(x - x2)*(x - x3)+A1*(x - x0)*(x - x2)*(x - x3)+
//A2*(x - x0)*(x - x1)*(x - x3)+A3*(x - x0)*(x - x1)*(x - x2)
let expr = some |> Infix.parseOrUndefined
let exp = expr |> Algebraic.expand
exp
|> Infix.format
|> printfn "%s"
(*
A0*x^3 + A1*x^3 + A2*x^3 + A3*x^3 - A1*x^2*x0 - A2*x^2*x0 - A3*x^2*x0 - A0*x^2*x1 -
A2*x^2*x1 - A3*x^2*x1 + A2*x*x0*x1 + A3*x*x0*x1 - A0*x^2*x2 - A1*x^2*x2 - A3*x^2*x2 +
A1*x*x0*x2 + A3*x*x0*x2 + A0*x*x1*x2 + A3*x*x1*x2 - A3*x0*x1*x2 - A0*x^2*x3 - A1*x^2*x3 -
A2*x^2*x3 + A1*x*x0*x3 + A2*x*x0*x3 + A0*x*x1*x3 + A2*x*x1*x3 - A2*x0*x1*x3 + A0*x*x2*x3 +
A1*x*x2*x3 - A1*x0*x2*x3 - A0*x1*x2*x3
*)
Раз уж речь зашла о полиномах Лагранжа приведу пример его вычисления по заданным
узлам интерполяции (x, y).
В качестве тестовых значений возьмем соответствующие из статьи в вики посвященной
интерполяционному многочлену Лагранжа.
Функция для создания "шаблонного выражения":
let generate n =
let y = sprintf "y%i"
let l i n = seq {
for j in 0..n - 1 do
if j <> i then
yield sprintf "(x - x%i)/(x%i - x%i)" j i j
}
let p i =
l i n
|> String.concat "*"
|> sprintf "%s*%s" <| y i
Array.init n p
|> String.concat "+"
задаем значения в виде словаря
let maps =
[
"x0", -1.5; "x1", -0.75; "x2", 0.0; "x3",0.75; "x4",1.5
"y0", -14.1014; "y1",-0.931596;"y2", 0.0;"y3",0.931596; "y4",14.1014
]
|> Map.ofList
Заменяем переменные на значения пробегая по словарю
let poly =
maps
|> Map.fold
(fun (acc : string) key value -> acc.Replace(key, string value))
(generate n)
Преобразуем в выражение и переведем в более читабельный вид:
let ex = poly |> Infix.parseOrUndefined
ex
|> Algebraic.expand
|> Infix.format
|> printfn "%s"
В результате получим следующее выражение:
(-1.47747377777778)*x + 4.83484760493827*x^3
Если у вас возникнут какие-либо вопросы связанные с F# или с библиотекой MathNet.Symbolics
не стесняйтесь спрашивать. Для удобства можете пинговать меня в F# чате на SO.
Ответ 2
Python + sympy
import sympy
x = sympy.Symbol('x')
xn = sympy.symarray('x', 3) # от x_0 до x_2
an = sympy.symarray('a', 3) # от a_0 до a_2
F = an[0]*(x-xn[0])*(x-xn[1])*(x-xn[2])+an[1]*(x-xn[0])*(x-xn[1])*(x-xn[2])+an[2]*(x-xn[0])*(x-xn[1])*(x-xn[2])
print(F.as_poly())
Результат:
Poly(x**3*a_0 + x**3*a_1 + x**3*a_2 - x**2*a_0*x_0 - x**2*a_0*x_1 - x**2*a_0*x_2
- x**2*a_1*x_0 - x**2*a_1*x_1 - x**2*a_1*x_2 - x**2*a_2*x_0 - x**2*a_2*x_1 - x**2*a_2*x_2
+ x*a_0*x_0*x_1 + x*a_0*x_0*x_2 + x*a_0*x_1*x_2 + x*a_1*x_0*x_1 + x*a_1*x_0*x_2 + x*a_1*x_1*x_2
+ x*a_2*x_0*x_1 + x*a_2*x_0*x_2 + x*a_2*x_1*x_2 - a_0*x_0*x_1*x_2 - a_1*x_0*x_1*x_2
- a_2*x_0*x_1*x_2, x, a_0, a_1, a_2, x_0, x_1, x_2, domain='ZZ')
Ответ 3
Рассмотрим (на примере PHP), как можно добиться результата без подключения внешних
библиотек.
Полином Лагранжа представляет собой сумму частичных полиномов Lk(x), каждый из которых
обращается в нуль во всех узлах, кроме узла с индексом k. Коэффициенты ak при частичных
полиномах принято записывать в виде отношения yk / Pk(xk), причём yk - значение полинома
в этом узле.
Коэффициентами каждого из частичных полиномов Лагранжа являются значения симметрических
полиномов в узлах, для вычисления которых можно написать короткую рекурсивную функцию
get_symm($ar, $k).
В дальнейшей работе можно опираться на функцию poly($ar, $factor), в которой обработка
одномерного массива $ar зависит от типа параметра $factor:
Если это число с плавающей точкой, то массив $ar на него умножается.
Если это такой же одномерный массив, то элементы этих массивов перемножаются.
Если $factor - это двумерный массив с той же "внешней" размерностью, что и массив
$ar, то вычисляется линейная комбинация его одномерных массивов, коэффициенты которой
берутся из массива $ar.
Если $factor - это целое число, то вычисляются коэффициенты частичного полинома Лагранжа
с таким индексом.
А если $factor - строка, то она задаёт функцию над массивом $ar.
Реализованы следующие функции:
'symmetric'. Вычисление массива значений всех возможных симметрических полиномов,
образуемых числами массива $ar.
'denominators'. Вычисление массива из всех знаменателей Pk(xk).
'degrees'. Вычисление квадратной матрицы, образованной степенями точек массива.
'reduced'. Вычисление коэффициентов полинома по его корням.
Пользовательская функция над массивом $ar, заданная вне функции poly().
Перечень опций с примерами выводится в начале работы программы. Такой подход, помимо
решения и тестирования поставленной задачи, даёт удобный инструмент для расширения
функциональности (интегрирование и дифференцирование полиномов, представление данных
и т.д.).
Программа на PHP:
$points = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0];
$coeffs = [-8.68, 83.44, -291.93, 398.63, -196.79, 26.04];
function print_m($text, $arr, $level=0){
$space = str_repeat(" ", $level++);
echo "$space$text";
if(gettype($arr)!="array"){
var_dump($arr);
return;
}
$flag = false;
foreach($arr as $value) $flag = $flag || (gettype($value)=="array");
foreach($arr as $key => $value) {
if(gettype($value) != "array") {
echo $flag ? "
$key => $value" : " $value";
} else {
print_m("
$space$key => [", $value, $level);
echo " ]";
}
}
$level--;
}
function get_symm($arr, $k){ // $k-th symmetric polynomial of $arr
if($k==1) return array_sum($arr);
$el = array_shift($arr); // Shifting of the first value from array
$cnt = count($arr);
return ($cnt == 1) ? ($el*end($arr)):
($cnt == $k-1) ? ($el*get_symm($arr, $k-1)):
$el * get_symm($arr, $k-1) + (get_symm($arr, $k));
}
function squares($a){
return $a * $a;
}
function round8($a){
return round($a, 8);
}
function poly($ar, $factor = ""){
switch (gettype($factor)) {
case 'float':
case 'double': // factoring
return array_map(function($val)use($factor){ return $val * $factor;},
$ar);
case 'array':
if(!is_array(end($factor))){ // arrays factoring
return array_map(function($a, $b){return $a * $b;}, $ar, $factor);
} else { // weighted sum of $factor
$sum = array_fill(0, count($ar), 0.0); // sum initialization
array_map(function($a, $ff)use(&$sum){
$f = poly($ff, $a); // array factoring
array_walk($f, function($val, $key)use(&$sum){$sum[$key] += $val;});
// summation
}, $ar, $factor);
return $sum;
}
case 'integer': // k-th Lagrange polynomial by roots
unset($ar[$factor]);
return poly($ar, 'reduced');
case 'string':
switch ($factor) {
case 'symmetric': // all symmetric polynomials's values
return array_map(function($a)use($ar){return get_symm($ar, $a);},
range(1,count($ar)));
case 'denominators':
$prod = []; // prod initialization
array_walk($ar, function($value, $key)use($ar, &$prod){
$p = 1;
array_walk($ar, function($val, $k)use($value, $key, &$p){
if($k != $key) $p *= $value - $val;
});
$prod[] = $p;
});
return $prod;
case 'degrees': // degrees in accordance with the point quantity
$cnt = count($ar);
$degree = array_fill(0, $cnt, 1.0);
$degrees = [$degree];
while(count($degrees) < $cnt) $degrees[] = ($degree = poly($degree,
$ar));
return $degrees;
case 'reduced': // all reduced polynomial's coefficients by roots
return array_merge(array_reverse(poly(poly($ar, -1.0),"symmetric")),[1.0]);
default:
return array_map($factor, $ar);
}
default:
return $arr;
}
}
$test = [1.0, 2.0, 3.0, 4.0];
$results[] = ["Testing points (in float point form) are: \$test " => $test];
$results[] = ["Array's factoring: poly(\$test, 5.0) " => poly($test, 5.0)];
$results[] = ["Multiplying to 1d array: poly(\$test, array_reverse(\$test))
" => poly($test, array_reverse($test))];
$results[] = ["Symmetric polynomials calculation: poly(\$test, 'symmetric')
" => poly($test, 'symmetric')];
$results[] = ["Lagrange denominators calculation: poly(\$test, 'denominators')
" => poly($test, 'denominators')];
$results[] = ["Square matrix of degrees: \$degrees = poly(\$test, 'degrees')
" => $degrees = poly($test, 'degrees')];
$results[] = ["2d array weighting via scalar production: poly(\$test, \$degrees)"
=> poly($test, $degrees)];
$results[] = ["Reduced polynomials via its roots: poly(\$test, 'reduced') "
=> poly($test, 'reduced')];
$results[] = ["User function appying: poly(\$test, 'squares') " => poly($test,
'squares')];
foreach ($test as $key => $value) {
$results[] = ["Lagrange partial polynomials: poly(\$test, $key) " => poly($test,
$key)];
}
print_m("Poly Functions Description:
", $results);
$symp[] = ["Lagrange Polynomial's Nodes " => $points];
$symp[] = ["Lagrange Partial Polynomials' Coefficients " => $coeffs];
$parts = [];
foreach ($points as $key =>$value){
$part = poly($points, $key);
$parts[] = $part;
$symp[] = ["Lagrange Partial Polynomial \#$key " => $part];
}
$lagrange = poly($coeffs, $parts);
$symp[] = ["Lagrange Polynomial's Coefficients " => poly($lagrange, 'round8')];
print_m("
Lagrange Polynomial Symplifying:
", $symp);
$denoms = poly($points, 'denominators');
$checking[] = ["Lagrange Denominators " => $denoms];
$values = poly(poly($coeffs, $denoms), 'round8');
$checking[] = ["Polynomial's Values In The Nodes Via Denominators" => $values];
$degrees = poly($points, 'degrees');
$checking[] = ["Degrees of issue Points " => $degrees];
$values2 = poly(poly($lagrange, $degrees), 'round8');
$checking[] = ["Polynomial's Values In The Nodes Via Symplified Polinomial" => $values2];
print_m("
Checking:
", $checking);
Результаты:
Poly Functions Description:
0 => [
Testing points (in float point form) are: $test => [ 1 2 3 4 ] ]
1 => [
Array's factoring: poly($test, 5.0) => [ 5 10 15 20 ] ]
2 => [
Multiplying to 1d array: poly($test, array_reverse($test)) => [ 4 6 6 4 ] ]
3 => [
Symmetric polynomials calculation: poly($test, 'symmetric') => [ 10 35 50 24 ] ]
4 => [
Lagrange denominators calculation: poly($test, 'denominators') => [ -6 2 -2 6 ] ]
5 => [
Square matrix of degrees: $degrees = poly($test, 'degrees')
=> [
0 => [ 1 1 1 1 ]
1 => [ 1 2 3 4 ]
2 => [ 1 4 9 16 ]
3 => [ 1 8 27 64 ] ] ]
6 => [
2d array weighting via scalar production: poly($test, $degrees) => [ 10 49 142 313 ] ]
7 => [
Reduced polynomials via its roots: poly($test, 'reduced') => [ 24 -50 35 -10 1 ] ]
8 => [
User function appying: poly($test, 'squares') => [ 1 4 9 16 ] ]
9 => [
Lagrange partial polynomials: poly($test, 0) => [ -24 26 -9 1 ] ]
10 => [
Lagrange partial polynomials: poly($test, 1) => [ -12 19 -8 1 ] ]
11 => [
Lagrange partial polynomials: poly($test, 2) => [ -8 14 -7 1 ] ]
12 => [
Lagrange partial polynomials: poly($test, 3) => [ -6 11 -6 1 ] ]
Lagrange Polynomial Symplifying:
0 => [
Lagrange Polynomial's Nodes => [ 0 0.2 0.4 0.6 0.8 1 ] ]
1 => [
Lagrange Partial Polynomials' Coefficients => [ -8.68 83.44 -291.93 398.63 -196.79 26.04 ] ]
2 => [
Lagrange Partial Polynomial \#0 => [ -0.0384 0.4384 -1.8 3.4 -3 1 ] ]
3 => [
Lagrange Partial Polynomial \#1 => [ -0 0.192 -1.232 2.84 -2.8 1 ] ]
4 => [
Lagrange Partial Polynomial \#2 => [ -0 0.096 -0.856 2.36 -2.6 1 ] ]
5 => [
Lagrange Partial Polynomial \#3 => [ -0 0.064 -0.624 1.96 -2.4 1 ] ]
6 => [
Lagrange Partial Polynomial \#4 => [ -0 0.048 -0.488 1.64 -2.2 1 ] ]
7 => [
Lagrange Partial Polynomial \#5 => [ -0 0.0384 -0.4 1.4 -2 1 ] ]
8 => [
Lagrange Polynomial's Coefficients => [ 0.333312 1.256224 -0.4096 13.538 -24.428 10.71 ] ]
Checking:
0 => [
Lagrange Denominators => [ -0.0384 0.00768 -0.00384 0.00384 -0.00768 0.0384 ] ]
1 => [
Polynomial's Values In The Nodes Via Denominators => [ 0.333312 0.6408192 1.1210112 1.5307392 1.5113472 0.999936 ] ]
2 => [
Degrees of issue Points => [
0 => [ 1 1 1 1 1 1 ]
1 => [ 0 0.2 0.4 0.6 0.8 1 ]
2 => [ 0 0.04 0.16 0.36 0.64 1 ]
3 => [ 0 0.008 0.064 0.216 0.512 1 ]
4 => [ 0 0.0016 0.0256 0.1296 0.4096 1 ]
5 => [ 0 0.00032 0.01024 0.07776 0.32768 1 ] ] ]
3 => [
Polynomial's Values In The Nodes Via Symplified Polinomial => [ 0.333312 0.6408192 1.1210112 1.5307392 1.5113472 0.999936 ] ]
Ответ 4
Спасибо всем за помощь. Особенно @FoggyFinder и @insolor.
Выкладываю завершённый и рабочий вариант.
Программа по узлам и заданной функции (double f(double x){}) составляет многочлен
Лагранжа, приводит его к каноническому виду и вычисляет определённый интеграл (работает,
только для числовых значений).
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using MathNet.Symbolics;
namespace Polynom
{
class Program
{
static void Main(string[] args)
{
double[] x = { 0, 0.2, 0.4, 0.6, 0.8, 1.0 }; // n - определяется как
количество узлов минус один.
//double[] x = { -1.5, -0.75, 0, 0.75, 1.5 };
double[] coef = GetCoefLagrange(x);
string exp = GetPolynom(coef, x);
string simpExp = SimplifyPolynom(exp);
var defInt = GetDefInt(GetSimplifyCoef(simpExp), 0, 1);
Console.WriteLine("Многочлен Лагранжа:");
Console.WriteLine(exp);
Console.WriteLine("---------------------");
Console.WriteLine("Канонический вид:");
Console.WriteLine(simpExp);
Console.WriteLine("---------------------");
Console.WriteLine("Определённый интеграл от 0 до 1:");
Console.WriteLine(defInt);
Console.ReadLine();
}
static double f(double x)
{
//return x * x * x * x;
//return Math.Tan(x);
return Math.Exp(Math.Pow(x, 1 / 3) * Math.Sin(Math.PI * x)) / (2 + Math.Cos(Math.PI
* x));
}
static double[] GetCoefLagrange(double[] x)
{
int n = x.Length - 1;
double[] a = new double[n + 1];
for (int i = 0; i < n + 1; i++)
{
double den = 1;
for (int j = 0; j < n + 1; j++)
{
if (i != j)
{
den *= (x[i] - x[j]);
}
}
a[i] = f(x[i]) / den;
}
return a;
}
static string GetPolynom(double[] a, double[] x)
{
string exp = "";
int n = x.Length - 1;
var culture = System.Globalization.CultureInfo.InvariantCulture;
for (int i = 0; i < n + 1; i++)
{
if (a[i] < 0) { exp += "("; }
exp += a[i].ToString(culture);
if (a[i] < 0) { exp += ")"; }
for (int j = 0; j < n + 1; j++)
{
if (j != i)
{
exp += "*(x - ";
if (x[j] < 0) { exp += "("; }
exp += x[j].ToString(culture);
if (x[j] < 0) { exp += ")"; }
exp += ")";
}
}
if (i != n) { exp += " + "; }
}
return exp;
}
static string SimplifyPolynom(string exp)
{
return Infix.Format(Algebraic.Expand(Infix.ParseOrUndefined(exp)));
}
static double[] GetSimplifyCoef(string exp)
{
string[] membr = exp.Split(" + ".ToCharArray(), StringSplitOptions.RemoveEmptyEntries);
double[] coef = new double[membr.Count() + 1];
for (int i = 0; i < membr.Count(); i++)
{
string str = membr[i].Split("*".ToCharArray(), StringSplitOptions.RemoveEmptyEntries)[0].Replace('.',
',');
string xPow = (membr[i].Split("*".ToCharArray(), StringSplitOptions.RemoveEmptyEntries).Length
> 1)? membr[i].Split("*".ToCharArray(), StringSplitOptions.RemoveEmptyEntries)[1] : "";
int pow = 0;
if (!xPow.Contains('x')) { pow = 0; }
if (xPow == "x") { pow = 1; }
if (xPow.Contains('x') && xPow != "x") { pow = int.Parse(xPow[xPow.Length
- 1].ToString()); }
if (str.Contains("("))
{
str = str.Substring(1, str.Length - 2);
coef[pow] = Math.Round(double.Parse(str), 2);
}
else
{
coef[pow] = Math.Round(double.Parse(str), 2);
}
pow++;
}
return coef;
}
static double GetDefInt(double[] A, double from, double to)
{
string str = "";
var culture = System.Globalization.CultureInfo.InvariantCulture;
for (int i = 0; i < A.Count(); i++)
{
str += (A[i] / (i + 1)).ToString(culture) + "*x^" + (i + 1);
if (i != A.Count() - 1) { str += " + "; }
}
var exp = Infix.ParseOrUndefined(str);
var x = new Dictionary { { "x", 0.0 } };
var resFrom = Evaluate.Evaluate(x, exp);
x = new Dictionary { { "x", 1.0 } };
var resTo = Evaluate.Evaluate(x, exp).RealValue;
return resTo - resFrom.RealValue;
}
}
}
Вывод консоли:
Многочлен Лагранжа:
(-8.68055555555555)*(x - 0.2)*(x - 0.4)*(x - 0.6)*(x - 0.8)*(x - 1) + 83.4365435984129*(x
- 0)*(x - 0.4)*(x - 0.6)*(x - 0.8)*(x - 1) + (-291.931019062828)*(x - 0)*(x - 0.2)*(x
- 0.6)*(x - 0.8)*(x - 1) + 398.628301975219*(x - 0)*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x
- 1) + (-196.79094312252)*(x - 0)*(x - 0.2)*(x - 0.4)*(x - 0.6)*(x - 1) + 26.0416666666667*(x
- 0)*(x - 0.2)*(x - 0.4)*(x - 0.6)*(x - 0.8)
---------------------
Канонический вид:
0.333333333333333 + 1.2551290418413*x + (-0.40261625087755)*x^2 + 13.5213484261595*x^3
+ (-24.4111890498517)*x^4 + 10.7039944993951*x^5
---------------------
Определённый интеграл от 0 до 1:
1,108
