Пытаясь заставить Mathematica приблизить интеграл

Question

Я пытаюсь заставить Mathematica аппроксимировать интеграл, который является функцией различных параметров.Мне не нужно, чтобы он был очень точным - ответ будет дробным, и 5 цифр было бы неплохо, но я бы согласился всего на 2.
Проблема в том, что похоронен символический интегралв основном интеграле, и я не могу использовать NIntegrate на нем, так как его символическое значение.

    F[x_, c_] := (1 - (1 - x)^c)^c;
    a[n_, c_, x_] := F[a[n - 1, c, x], c];
    a[0, c_, x_] = x;

    MyIntegral[n_,c_] := 
      NIntegrate[Integrate[(D[a[n,c,y],y]*y)/(1-a[n,c,x]),{y,x,1}],{x,0,1}]

Mathematica начинает зависать, когда n больше 2 и c больше 3 или около тогокак и n, и c немного выше).

Существуют ли какие-либо приемы для переписывания этого выражения, чтобы его было легче оценить?Я играл с различными опциями WorkingPrecision и AccuracyGoal и PrecisionGoal на внешней NIntegrate, но ничего из этого не помогает внутреннему интегралу, в котором проблема.На самом деле, для более высоких значений n и c я даже не могу заставить Mathematica расширить внутреннюю производную, то есть зависает

Expand[D[a[4,6,y],y]] 

.

Я использую Mathematica 8 для студентов.

Если у кого-нибудь есть какие-либо советы о том, как я могу заставить М. приблизиться к этому, я был бы признателен.

Answer 1

Поскольку вам нужен только числовой вывод (или это то, что вы получите в любом случае), вы можете преобразовать символьную интеграцию в числовую, используя просто NIntegrate следующим образом:

Clear[a,myIntegral]
a[n_Integer?Positive, c_Integer?Positive, x_] := 
  a[n, c, x] = (1 - (1 - a[n - 1, c, x])^c)^c;
a[0, c_Integer, x_] = x;

myIntegral[n_, c_] := 
 NIntegrate[D[a[n, c, y], y]*y/(1 - a[n, c, x]), {x, 0, 1}, {y, x, 1},
   WorkingPrecision -> 200, PrecisionGoal -> 5]

Этогораздо быстрее, чем выполнять символическую интеграцию.Вот сравнение:

йода:

myIntegral[2,2]//Timing
Out[1]= {0.088441, 0.647376595...}

myIntegral[5,2]//Timing
Out[2]= {1.10486, 0.587502888...}

rcollyer:

MyIntegral[2,2]//Timing
Out[3]= {1.0029, 0.647376}

MyIntegral[5,2]//Timing 
Out[4]= {27.1697, 0.587503006...}
(* Obtained with WorkingPrecision->500, PrecisionGoal->5, MaxRecursion->20 *)

Функция Jand имеет тайминги, подобные rcollyer.Конечно, когда вы увеличиваете n, вам придется увеличивать свой WorkingPrecision намного выше, чем , который вы испытали в своем предыдущем вопросе .Поскольку вы сказали, что вам нужно только около 5 цифр точности, я явно установил PrecisionGoal на 5. Вы можете изменить это в соответствии со своими потребностями.

Answer 2

Чтобы систематизировать комментарии, я бы попробовал следующее. Во-первых, чтобы исключить бесконечную рекурсию в отношении переменной n, я бы переписал ваши функции как

F[x_, c_] := (1 - (1-x)^c)^c;
(* see note below *)
a[n_Integer?Positive, c_, x_] := F[a[n - 1, c, x], c];  
a[0, c_, x_] = x;

таким образом n==0 будет фактически точкой остановки. Форма ?Positive является PatternTest и полезна для применения дополнительных условий к параметрам. Я подозреваю, что проблема в том, что NIntegrate переоценивает внутреннюю Integrate для каждого значения x, поэтому я бы извлек эту оценку, как

MyIntegral[n_,c_] := 
  With[{ int = Integrate[(D[a[n,c,y],y]*y)/(1-a[n,c,x]),{y,x,1}] },
    NIntegrate[int,{x,0,1}]
  ]

где With - одна из нескольких областей видимости, специально предназначенных для создания локальных констант.

Ваши комментарии указывают на то, что внутренний интеграл занимает много времени, пытались ли вы упростить подынтегральное выражение, поскольку оно является производной от a от функции a? Мне кажется, это результат расширения правила цепочки.

Примечание : согласно предложению Йоды в комментариях, вы можете добавить механизм кэширования или запоминания к a. Измените его определение на

d:a[n_Integer?Positive, c_, x_] := d = F[a[n - 1, c, x], c];

Хитрость в том, что в d:a[ ... ], d является именованным шаблоном, который снова используется в d = F[...] для кэширования значения a для этих конкретных значений параметров.