Erreur d'arrondi Scilab

Question

Je ne peux pas résoudre un problème dans Scilab car il est bloqué à cause d'erreurs d'arrondi. Je reçois le message

!--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) 

C'est un processus Gram Schmidt avec l'intégrale du produit de deux fonctions et un poids comme le produit scalaire, entre -1 et 1. gram_schmidt_pol est le processus spécialement conçu pour polynome et scalpol est le produit scalaire décrit pour le polynôme.

Le a et b sont des paramètres pour le weigth, qui est (1+x)^a*(1-x)^b

L'entrée est une matrice qui représente un ensemble de vecteurs, il fonctionne bien avec la matrice [[1;2;3],[4;5;6],[7;8;9]], mais il échoue avec l'erreur de message ci-dessus sur la matrice eye(2,2), en plus de cela, je dois le faire sur les yeux (9,9)!

J'ai cherché un "paramètre de tolérance" dans les menus, il y en a dans General->Preferences->Xcos->Simulation mais je crois que ce n'est pas pour ce que je voulais, j'ai essayé les réglages bas (haute tolérance) dedans et ça n'a pas t changer quoi que ce soit.

Alors, comment puis-je résoudre ce problème de rebond?

N'hésitez pas à me dire que mon message manque de clarté. Merci.

Modifier: Code des fonctions:

// 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

Par défaut, la routine de Scilab integrate tente de réaliser une erreur absolue au plus 1E-8 et erreur relative au plus 1E-14. C'est raisonnable, mais son traitement de l'erreur relative ne prend pas en compte les problèmes qui se produisent lorsque la valeur exacte est zéro. (Voir How to calculate relative error when true value is zero?). Pour cette raison, même le simple

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

génère une erreur (dans Scilab 5.5.1).

Et c'est ce qui se passe dans le processus d'exécution de votre programme: certaines intégrales sont nulles. Il y a deux solutions:

(A) Renoncer à l'erreur relative liée, en spécifiant comme 1:

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

(B) Ajouter une constante à la fonction en cours d'intégration, puis soustrayez du résultat :

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

(Ce dernier devrait travailler dans la plupart des cas, mais si l'intégrale se trouve être exactement -200, vous aurez le même problème)

---

Ce qui précède fonctionne pour gram_schmidt_pol(eye(2,2), -1/2, -1/2) mais pour plus grand, par exemple, gram_schmidt_pol(eye(9,9), -1/2, -1/2), il renvoie l'erreur "L'intégrale est probablement divergente, ou lentement convergente".

Il semble que la routine d'intégration adaptative ne puisse pas gérer les fonctions du type que vous avez. Une solution de rechange consiste à utiliser le inttrap simple, qui applique simplement la règle trapézoïdale. Comme en x = -1 et 1, la fonction poids n'est pas définie, les extrémités doivent être exclues.

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 

Pour que cela fonctionne, d'autres fonctions connexes doivent être vectorisés (* et^changé pour * et^le cas échéant..):

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 

Le résultat est garanti pour fonctionner, bien que la précision peut être un peu plus faible: vous allez obtenir des nombres comme 3D-16 qui sont en réalité des zéros.