Я могу предложить способ свести ваше уравнение к интегральному уравнению, которое может быть решено численно путем приближения его ядра к матрице, тем самым сводя интеграцию к умножению матриц.
Во-первых, ясно, что уравнение может быть интегрировано дважды по x, сначала от 1 до x, а затем от 0 до x, так что:
Теперь мы можем дискретизировать это уравнение, поместив его на эквидистантную сетку:
Здесь A[x] становится вектороми интегрированное ядро iniIntK становится матрицей, а интеграция заменяется умножением матрицы.Затем задача сводится к системе линейных уравнений.
Самый простой случай (который я рассмотрю здесь) - это когда ядро iniIntK может быть получено аналитически - в этом случае этот метод будет довольно быстрым.Вот функция для создания интегрированного ядра в виде чистой функции:
Clear[computeDoubleIntK]
computeDoubleIntK[kernelF_] :=
Block[{x, x1},
Function[
Evaluate[
Integrate[
Integrate[kernelF[y, x1], {y, 1, x}] /. x -> y, {y, 0, x}] /.
{x -> #1, x1 -> #2}]]];
В нашем случае:
In[99]:= K[x_,x1_]:=1;
In[100]:= kernel = computeDoubleIntK[K]
Out[100]= -#1+#1^2/2&
Вот функция для получения матрицы ядра и rh, svector:
computeDiscreteKernelMatrixAndRHS[intkernel_, a0_, aprime1_ ,
delta_, interval : {_, _}] :=
Module[{grid, rhs, matrix},
grid = Range[Sequence @@ interval, delta];
rhs = a0 + aprime1*grid; (* constant plus a linear term *)
matrix =
IdentityMatrix[Length[grid]] - delta*Outer[intkernel, grid, grid];
{matrix, rhs}]
Чтобы дать очень грубое представление о том, как это может выглядеть (я использую здесь delta = 1/2):
In[101]:= computeDiscreteKernelMatrixAndRHS[kernel,0,1,1/2,{0,1}]
Out[101]= {{{1,0,0},{3/16,19/16,3/16},{1/4,1/4,5/4}},{0,1/2,1}}
Теперь нам нужно решить линейное уравнение, иинтерполировать результат, что делается с помощью следующей функции:
Clear[computeSolution];
computeSolution[intkernel_, a0_, aprime1_ , delta_, interval : {_, _}] :=
With[{grid = Range[Sequence @@ interval, delta]},
Interpolation@Transpose[{
grid,
LinearSolve @@
computeDiscreteKernelMatrixAndRHS[intkernel, a0, aprime1, delta,interval]
}]]
Здесь я назову его с delta = 0.1:
In[90]:= solA = computeSolution[kernel,0,1,0.1,{0,1}]
Out[90]= InterpolatingFunction[{{0.,1.}},<>]
Теперь мы построим график зависимости результата от точногоаналитическое решение, найденное @Sasha, а также ошибка:
Я специально выбрал delta достаточно большим, чтобы ошибки были видны.Если вы выбрали delta скажем 0.01, графики будут визуально идентичны.Конечно, цена взятия меньшего delta - это необходимость производить и решать большие матрицы.
Для ядер, которые могут быть получены аналитически, основное узкое место будет в LinearSolve, но на практике это довольно быстро (для матриц не слишком больших).Когда ядра не могут быть интегрированы аналитически, основным узким местом будет вычисление ядра во многих точках (создание матрицы. Обратная матрица имеет большую асимптотическую сложность, но это начнет играть роль для действительно больших матриц - которые не нужны вэтот подход, так как он может быть объединен с итеративным - см. ниже).Обычно вы определяете:
intK[x_?NumericQ, x1_?NumericQ] := NIntegrate[K[y, x1], {y, 1, x}]
intIntK[x_?NumericQ, x1_?NumericQ] := NIntegrate[intK[z, x1], {z, 0, x}]
В качестве способа ускорения в таких случаях вы можете предварительно вычислить ядро intK в сетке, а затем интерполировать, и то же самое для intIntK.Это, однако, приведет к дополнительным ошибкам, которые вам придется оценить (учесть).
Сама сетка не должна быть равноудаленной (я просто использовал ее для простоты), но может (и, вероятно, должна) быть адаптивнойи вообще неравномерно.
В качестве окончательной иллюстрации рассмотрим уравнение с нетривиальным, но символически интегрируемым ядром:
In[146]:= sinkern = computeDoubleIntK[50*Sin[Pi/2*(#1-#2)]&]
Out[146]= (100 (2 Sin[1/2 \[Pi] (-#1+#2)]+Sin[(\[Pi] #2)/2]
(-2+\[Pi] #1)))/\[Pi]^2&
In[157]:= solSin = computeSolution[sinkern,0,1,0.01,{0,1}]
Out[157]= InterpolatingFunction[{{0.,1.}},<>]
Вот несколько проверок:
In[163]:= Chop[{solSin[0],solSin'[1]}]
Out[163]= {0,1.}
In[153]:=
diff[x_?NumericQ]:=
solSin''[x] - NIntegrate[50*Sin[Pi/2*(#1-#2)]&[x,x1]*solSin[x1],{x1,0,1}];
In[162]:= diff/@Range[0,1,0.1]
Out[162]= {-0.0675775,-0.0654974,-0.0632056,-0.0593575,-0.0540479,-0.0474074,
-0.0395995,-0.0308166,-0.0212749,-0.0112093,0.000369261}
В заключение хочу подчеркнуть, что для этого метода необходимо провести тщательный анализ ошибок, чего я не делал.
РЕДАКТИРОВАТЬ
Вы также можете использовать этот метод, чтобы получить начальное приближенное решение, а затем итеративно улучшить его, используя FixedPoint или другими способами - таким образом, вы будете иметь относительно быструю сходимость и сможете достичь требуемой точности безнеобходимость строить и решать огромные матрицы.