wolfram-mathematica - 如何在 NDSolve 中引用函数的特定点？

Question

问题:
我正在尝试解决这个微分方程:

K[x_, x1_] := 1;
NDSolve[{A''[x] == Integrate[K[x, x1] A[x1], {x1, 0, 1}], 
         A[0] == 0, A'[1] == 1}, A[x], x]

我收到错误( Function::slotn 和 NDSolve::ndnum )
(它应该返回一个等于 3/16 x^2 + 5/8 x 的数字函数)
我正在寻找一种方法来解决这个微分方程:有没有办法以更好的形式写出来，这样 NDSolve 就会理解它？是否有其他功能或包可以提供帮助？
注 1: 在我的完整问题中， K[x, x1] 不是 1——它(以复杂的方式)取决于 x 和 x1 。
注 2: 天真地推导出方程的两侧关于 x 是行不通的，因为积分限制是确定的。
我的第一印象:
似乎 Mathematica 不喜欢我在 A[x] 中引用一个点——当我在做这个简化版本时会发生同样的错误:

NDSolve[{A''[x] == A[0.5], A[0] == 0, A'[1] == 1}, A[x], x]

(它应该返回一个等于 2/11 x^2 + 7/11 x 的数字函数)
在这种情况下，可以通过解析求解 A''[x] == c ，然后找到 c 来避免这个问题，但在我的第一个问题中它似乎不起作用——它只是将微分方程转换为一个积分方程，而 (N)DSolve 不能解决然后。

Answer 1

我可以建议一种将方程简化为积分方程的方法，该方程可以通过将其内核与矩阵近似来进行数值求解，从而减少对矩阵乘法的积分。

首先，显然方程可以在 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&

这是产生核矩阵和 r.h,s 向量的函数:

 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或者其他方式迭代改进——这样你就会有比较快的收敛，不需要构造求解就能达到要求的精度巨大的矩阵。