Ошибка округления Scilab

Question

Я не могу решить проблему в Scilab, потому что она застревает из-за ошибок округления. Я получаю сообщение

!--error 9999 
Error: Round-off error detected, the requested tolerance (or default) cannot be achieved. Try using bigger tolerances. 
at line  2 of function scalpol called by : 
at line  7 of function gram_schmidt_pol called by : 
gram_schmidt_pol(a,-1/2,-1/2) 

Это процесс Грама Шмидта с интегралом от произведения двух функций и весом, как скалярное произведение между -1 и 1. gram_schmidt_pol процесса разработан специально для полинома, и scalpol - скалярный продукт, описанный для полинома.

a и b являются параметрами для Weigth, который (1+x)^a*(1-x)^b

Запись является матрица, представляющая собой набор векторов, он хорошо работает с матрицей [[1;2;3],[4;5;6],[7;8;9]], но он терпит неудачу с вышеуказанной ошибкой сообщения на матрице eye(2,2), в дополнение к этому, мне нужно сделать это на глаз (9,9)!

Я искал «настройку допуска» в меню, есть некоторые в General->Preferences->Xcos->Simulation, но я считаю, что это не то, что я не хочу, я пробовал низкие настройки (высокая толерантность) в нем, и он hasn ' ничего не меняйте.

Так как я могу решить эту проблему с rounf-off?

Не стесняйтесь сказать мне, что мое сообщение не имеет ясности. Спасибо.

Редактировать: Код функции:

// function that evaluate a polynomial (vector of coefficients) in x 
function [y] = pol(p, x) 
    y = 0 
    for i=1:length(p) 
     y = y + p(i)*x^(i-1) 
    end 
endfunction 

// weight function evaluated in x, parametrized by a and b 
// (poids = weight in french) 
function [y] = poids(x, a, b) 
    y = (1-x)^a*(1+x)^b 
endfunction 

// scalpol compute scalar product between polynomial p1 and p2 
// using integrate, the weight and the pol functions. 
function [s] = scalpol(p1, p2, a, b) 
    s = integrate('poids(x,a, b)*pol(p1,x)*pol(p2,x)', 'x', -1, 1) 
endfunction 

// norm associated to scalpol 
function [y] = normscalpol(f, a, b) 
    y = sqrt(scalpol(f, f, a, b)) 
endfunction 

// finally the gram schmidt process on a family of polynome 
// represented by a matrix 
function [o] = gram_schmidt_pol(m, a, b) 
    [n,p] = size(m) 
    o(1:n) = m(1:n,1)/(normscalpol(m(1:n,1), a, b)) 
    for k = 2:p 
     s =0 
     for i = 1:(k-1) 
      s = s + (scalpol(o(1:n,i), m(1:n,k), a, b)/scalpol(o(1:n,i),o(1:n,i), a, b) .* o(1:n,i)) 
     end 
     o(1:n,k) = m(1:n,k) - s 
     o(1:n,k) = o(1:n,k) ./ normscalpol(o(1:n,k), a, b) 
    end 
endfunction 

Answer 1

По умолчанию integrate процедура пытается Scilab для достижения абсолютной ошибки в большинстве 1e-8 и относительной погрешностью не более 1e-14. Это разумно, но его обработка относительной ошибки не учитывает проблемы, возникающие, когда точное значение равно нулю. (См. How to calculate relative error when true value is zero?). По этой причине даже простые

integrate('x', 'x', -1, 1) 

вызывает ошибку (в Scilab 5.5.1).

И это то, что происходит при запуске вашей программы: некоторые интегралы равны нулю. Есть два решения:

(А) Откажитесь от относительной погрешности, связанной, указав его в качестве 1:

integrate('...', 'x', -1, 1, 1e-8, 1) 

(B) Добавить некоторые константы к функциям интегрирован, а затем вычесть из результата :

integrate('100 + ... ', 'x', -1, 1) - 200 

(последний должен работать в большинстве случаев, хотя, если интеграл случается именно -200, вы будете иметь ту же проблему снова)

---

Вышеуказанные работы для gram_schmidt_pol(eye(2,2), -1/2, -1/2), но для большего, скажем, gram_schmidt_pol(eye(9,9), -1/2, -1/2), он выдает ошибку «Интеграл, вероятно, расходящийся или медленно сходящийся».

Похоже, что процедура адаптивной интеграции не может обрабатывать функции, которые у вас есть. Резервный должен использовать простой inttrap вместо этого, который просто применяет правило трапеции. Поскольку в x = -1 и 1 функция poids не определена, конечные точки должны быть исключены.

function [s] = scalpol(p1, p2, a, b) 
    t = -0.9995:0.001:0.9995 
    y = poids(t,a, b).*pol(p1,t).*pol(p2,t) 
    s = inttrap(t,y) 
endfunction 

Для того, чтобы это работало, другие связанные функции должны быть векторизации (* и^меняется на * и ^, где это необходимо..):

function [y] = pol(p, x) 
    y = 0 
    for i=1:length(p) 
     y = y + p(i)*x.^(i-1) 
    end 
endfunction 

function [y] = poids(x, a, b) 
    y = (1-x).^a.*(1+x).^b 
endfunction 

Результат гарантированно работать, хотя точность может быть немного ниже: вы получите некоторые цифры, такие как 3D-16, которые на самом деле являются нулями.