如何在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 - 它取决于(以复杂的方式)on 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&

这是生成内核矩阵和rh,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或以其他方式迭代地改进它- 这样您将获得相对快速的收敛并且能够达到所需的精度而无需构造和解决巨大的问题矩阵.