2014-09-18 4 views
1

Я пытаюсь интегрировать Спектр деления Ватта с помощью Sage (команды ниже), но когда я это сделал, я заметил, что результат интеграла не имеет смысла. Спектр деления Ватта - это pdf, а при интегрировании от 0 до 10 должен быть примерно один (подтвержденный Wolfram Alpha). С Sage, когда нижняя граница близка к 0, интеграл отрицателен, а затем прыгает положительно, когда нижняя граница ближе к 1.Неисправность шалфея на спектре деления Ватта

Ниже приведен график интеграла от спектра деления Ватта с границами x to 10 с x от 0 to 1. Скачок очень ясен.

Как исправить это и быть уверенным в будущем, что интегралы оцениваются правильно?

var('o') 
assume(x>10) 
plot(integral(0.453*exp(-1.036*o)*sinh(sqrt(2.29*o)), (o, x, 10)), (x,0,1)) 

Sage Plot

предложение Per Роберта Dodier, я попробовал интеграцию следующим образом:

map(var, 'abcde') 
assume(e-d>0) 
foo(a,b,c,d,e) = integral(a*exp(b*o)*sinh(sqrt(c*o)), (o, d, e)) 
float(foo(*map(QQ, [0.453,-1.036,2.29,0,10]))) 

Это еще возвращает значение -0.0010615612261540579.

ответ

0

Я не знаю, почему Sage не ретрансляции вопросов Maxima, но без них, вы не можете сказать (в Sage), что результат отличается в зависимости от параметров.

Вот вывод в Maxima. Я перефразировал интеграл как от 0 до x вместо x до 10. Так как подынтегральное выражение является pdf, мы должны видеть, что интеграл переходит от 0 до 1, так как увеличивается x.

Оригинальная формула.

(%i2) foo : 0.453*exp(-1.036*o)*sinh(sqrt(2.29*o)) $ 

Экспресс, что символы вместо цифр.

(%i3) foo1 : a * exp(-b*o) * sinh(sqrt(c*o)) $ 

Forestall несколько вопросов из integrate.Это не обязательно нужно делать; integrate просто задаст дополнительные вопросы, если эти предположения не указаны.

(%i4) assume (equal (a, ratsimp (0.453)), equal (b, ratsimp (1.036)), equal (c, ratsimp (2.29))) $ 

rat: replaced 0.453 by 453/1000 = 0.453 

rat: replaced 1.036 by 259/250 = 1.036 

rat: replaced 2.29 by 229/100 = 2.29 
(%i5) assume (x > 0) $ 

Вычислить интеграл. integrate спрашивает о x по отношению к параметрам a, b и c. Я дам каждый ответ по очереди и посмотрю, что я получу. I1, I2 и I3 - это результаты для всех ответов. Сначала я отвечу p (положительный).

(%i6) I1 : integrate (foo1, o, 0, x); 
Is 4*b^2*x-c positive, negative or zero? 

p; 
(%o6) a*(-gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
     *sqrt(c)*%e^(c/(4*b))*sqrt(x) 
     /(2*sqrt(b)*abs(2*b*sqrt(x)-sqrt(c))) 
     +gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
     *c*%e^(c/(4*b)) 
     /(4*b^(3/2)*abs(2*b*sqrt(x)-sqrt(c))) 
     -gamma_incomplete(1/2,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b)) 
     *sqrt(c)*%e^(c/(4*b)) 
     /(4*b^(3/2))+sqrt(%pi)*sqrt(c)*%e^(c/(4*b))/(2*b^(3/2)) 
     +gamma_incomplete(1,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b)) 
     *%e^(c/(4*b)) 
     /(2*b) 
     -gamma_incomplete(1,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
     *%e^(c/(4*b)) 
     /(2*b)) 

Теперь я отвечу n (отрицательный).

(%i7) I2 : integrate (foo1, o, 0, x); 
Is 4*b^2*x-c positive, negative or zero? 

n; 
(%o7) a*(-gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
     *sqrt(c)*%e^(c/(4*b))*sqrt(x) 
     /(2*sqrt(b)*abs(2*b*sqrt(x)-sqrt(c))) 
     +gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
     *c*%e^(c/(4*b)) 
     /(4*b^(3/2)*abs(2*b*sqrt(x)-sqrt(c))) 
     -gamma_incomplete(1/2,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b)) 
     *sqrt(c)*%e^(c/(4*b)) 
     /(4*b^(3/2)) 
     +gamma_incomplete(1,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b)) 
     *%e^(c/(4*b)) 
     /(2*b) 
     -gamma_incomplete(1,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
     *%e^(c/(4*b)) 
     /(2*b)) 

И вот z (ноль).

(%i8) I3 : integrate (foo1, o, 0, x); 
Is 4*b^2*x-c positive, negative or zero? 

z; 
(%o8) a*(-gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
     *sqrt(c)*%e^(c/(4*b))*sqrt(x) 
     /(2*sqrt(b)*abs(2*b*sqrt(x)-sqrt(c))) 
     +gamma_incomplete(1/2,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
     *c*%e^(c/(4*b)) 
     /(4*b^(3/2)*abs(2*b*sqrt(x)-sqrt(c))) 
     -gamma_incomplete(1/2,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b)) 
     *sqrt(c)*%e^(c/(4*b)) 
     /(4*b^(3/2)) 
     +gamma_incomplete(1,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b)) 
     *%e^(c/(4*b)) 
     /(2*b) 
     -gamma_incomplete(1,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
     *%e^(c/(4*b)) 
     /(2*b)) 

Я не могу сказать, просто посмотрев, возможно, некоторые из этих результатов одинаковы. Давай проверим.

(%i9) is (I1 = I2); 
(%o9) false 
(%i10) is (I1 = I3); 
(%o10) false 
(%i11) is (I2 = I3); 
(%o11) true 

ОК, так что есть два отличных результатов: I1 (ответил p) и I2 (ответил n и z). Я помещу эти два в одно условное выражение, в котором тест - это вопрос, который задал integrate. Обратите внимание, что c/(4*b^2) составляет около 0,533, что представляется точкой, в которой график, показанный OP, имеет скачок.

(%i12) I : if x < c/(4*b^2) then ''I2 else ''I1; 
(%o12) if x < c/(4*b^2) 
      then a*(-gamma_incomplete(1/2, 
            -(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
        *sqrt(c)*%e^(c/(4*b))*sqrt(x) 
        /(2*sqrt(b)*abs(2*b*sqrt(x)-sqrt(c))) 
        +gamma_incomplete(1/2, 
            -(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
        *c*%e^(c/(4*b)) 
        /(4*b^(3/2)*abs(2*b*sqrt(x)-sqrt(c))) 
        -gamma_incomplete(1/2,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b)) 
        *sqrt(c)*%e^(c/(4*b)) 
        /(4*b^(3/2)) 
        +gamma_incomplete(1,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b)) 
        *%e^(c/(4*b)) 
        /(2*b) 
        -gamma_incomplete(1,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
        *%e^(c/(4*b)) 
        /(2*b)) 
      else a*(-gamma_incomplete(1/2, 
            -(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
        *sqrt(c)*%e^(c/(4*b))*sqrt(x) 
        /(2*sqrt(b)*abs(2*b*sqrt(x)-sqrt(c))) 
        +gamma_incomplete(1/2, 
            -(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
        *c*%e^(c/(4*b)) 
        /(4*b^(3/2)*abs(2*b*sqrt(x)-sqrt(c))) 
        -gamma_incomplete(1/2,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b)) 
        *sqrt(c)*%e^(c/(4*b)) 
        /(4*b^(3/2))+sqrt(%pi)*sqrt(c)*%e^(c/(4*b))/(2*b^(3/2)) 
        +gamma_incomplete(1,(4*b^2*x+4*b*sqrt(c)*sqrt(x)+c)/(4*b)) 
        *%e^(c/(4*b)) 
        /(2*b) 
        -gamma_incomplete(1,-(-4*b^2*x+4*b*sqrt(c)*sqrt(x)-c)/(4*b)) 
        *%e^(c/(4*b)) 
        /(2*b)) 

Теперь я заговорю это. plot2d оценивается с помощью a, b и c, временно привязанных к определенным значениям. Участок выглядит так, как ожидалось.

(%i13) plot2d (I, [x, 0, 10]), a=0.453, b=1.036, c=2.29; 

enter image description here

+0

Wow thanks! Интересно, как работают эти предположения. Немного страшно, как легко разбить интеграл, когда вы не знаете, как эти системы работают в бэкэнд. Я попытаюсь поднять вопрос с людьми Sage, что сообщения не передаются от Maxima –

0

Sage, по умолчанию, передает Maxima для вычисления интегралов. Максимум сам по себе задает некоторые вопросы о признаках пары переменных и дает разные результаты, учитывая разные ответы на эти вопросы. Я думаю, что вопросы потерялись в интерфейсе Sage/Maxima. Мудрец может также вызвать Sympy для вычисления интегралов, но я обнаружил, что Sympy не может решить этот интеграл.

Имейте в виду, что Maxima более удобен с точными числами (т. Е. Интегралами и рациональными), чем неточные (т. Е. Плавает). Мой совет - заменить поплавки рациональными или символами (позже подставляя числовые значения для символов).

Учитывая разные результаты, вы можете скомпоновать функцию, похожую на if <condition 1> then <result 1> else if <condition 2> then <result 2> else .... В настоящее время нет хорошего способа сделать это автоматически. (Я работал над экспериментальным пакетом с именем noninteractive, но я не могу рекомендовать его для общего использования, более надежным является создание if--then.)

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

(%i2) foo : 0.453*exp(-1.036*o)*sinh(sqrt(2.29*o)) $ 
(%i3) bar : ratsimp (foo); 

rat: replaced -1.036 by -259/250 = -1.036 

rat: replaced 1.513274595042156 by 35244471/23290202 = 1.513274595042156 

rat: replaced 0.453 by 453/1000 = 0.453 
(%o3) 453*sinh(35244471*sqrt(o)/23290202)*%e^-(259*o/250)/1000 
(%i4) I1 : integrate (bar, o, x, 10); 
Is 36386982230699133124*x-19408949001091265625 positive, negative or zero? 

p; 
Is x-344460873305900065615/36386982230699133124 positive, negative or zero? 

p; 
(%o4) 453*(22027794375*sqrt(10)*%e^(155271592008730125/280980557766016472) 
         *gamma_incomplete(1/2, 
             -(-36386982230699133124*x 
             +53150092471010944500*sqrt(x) 
             -19408949001091265625) 
             /35122569720752059000)*sqrt(x) 
      /(2*sqrt(259)*abs(6032162318*sqrt(x)-4405558875)) 
      -97044745005456328125*sqrt(10) 
           *%e^(155271592008730125/280980557766016472) 
           *gamma_incomplete(1/2, 
               -(-36386982230699133124*x 
                +53150092471010944500 
                *sqrt(x) 
                -19408949001091265625) 
                /35122569720752059000) 
      /(46580404*259^(3/2)*abs(6032162318*sqrt(x)-4405558875)) 
      -125*%e^(155271592008730125/280980557766016472) 
       *gamma_incomplete(1, 
           (36386982230699133124*x 
           +53150092471010944500*sqrt(x) 
           +19408949001091265625) 
           /35122569720752059000) 
      /259 
      +125*%e^(155271592008730125/280980557766016472) 
       *gamma_incomplete(1, 
           -(-36386982230699133124*x 
           +53150092471010944500*sqrt(x) 
           -19408949001091265625) 
           /35122569720752059000) 
      /259 
      +291127525*10^(3/2)*%e^(155271592008730125/280980557766016472) 
        *gamma_incomplete(1, 
             (106300184942021889*10^(5/2) 
             +76655754261616519373) 
             /7024513944150411800) 
      /(6032162318*sqrt(10)-4405558875) 
      -550694859375*%e^(155271592008730125/280980557766016472) 
         *gamma_incomplete(1, 
             (106300184942021889*10^(5/2) 
              +76655754261616519373) 
              /7024513944150411800) 
      /(259*(6032162318*sqrt(10)-4405558875)) 
      -291127525*10^(3/2)*%e^(155271592008730125/280980557766016472) 
        *gamma_incomplete(1, 
             -(106300184942021889*10^(5/2) 
             -76655754261616519373) 
             /7024513944150411800) 
      /(6032162318*sqrt(10)-4405558875) 
      +550694859375*%e^(155271592008730125/280980557766016472) 
         *gamma_incomplete(1, 
             -(106300184942021889*10^(5/2) 
              -76655754261616519373) 
              /7024513944150411800) 
      /(259*(6032162318*sqrt(10)-4405558875)) 
      +22027794375*sqrt(10)*%e^(155271592008730125/280980557766016472) 
         *gamma_incomplete(1/2, 
             (36386982230699133124*x 
             +53150092471010944500*sqrt(x) 
             +19408949001091265625) 
             /35122569720752059000) 
      /(46580404*259^(3/2)) 
      -110138971875*%e^(155271592008730125/280980557766016472) 
         *gamma_incomplete(1/2, 
             (106300184942021889*10^(5/2) 
              +76655754261616519373) 
              /7024513944150411800) 
      /((6032162318*sqrt(10)-4405558875)*sqrt(259)) 
      +97044745005456328125*sqrt(10) 
           *%e^(155271592008730125/280980557766016472) 
           *gamma_incomplete(1/2, 
               (106300184942021889*10^(5/2) 
                +76655754261616519373) 
                /7024513944150411800) 
      /(46580404*(6032162318*sqrt(10)-4405558875)*259^(3/2)) 
      -110138971875*%e^(155271592008730125/280980557766016472) 
         *gamma_incomplete(1/2, 
             -(106300184942021889*10^(5/2) 
              -76655754261616519373) 
              /7024513944150411800) 
      /((6032162318*sqrt(10)-4405558875)*sqrt(259)) 
      +97044745005456328125*sqrt(10) 
           *%e^(155271592008730125/280980557766016472) 
           *gamma_incomplete(1/2, 
               -(106300184942021889*10^(5/2) 
                -76655754261616519373) 
                /7024513944150411800) 
      /(46580404*(6032162318*sqrt(10)-4405558875)*259^(3/2))) 
/1000 
+0

Спасибо за предложение об использовании символов. Я попробовал и добавил некоторые вопросы на свои вопросы. Извините, но я не совсем понял, что вы пытались сказать о структурах if-then. –

+0

@chewsocks Я добавил еще один ответ, который работает через материал с символами. –

 Смежные вопросы

  • Нет связанных вопросов^_^