方程式を積分方程式に縮小する方法を提案できます。これは、そのカーネルを行列で近似することで数値的に解くことができ、それによって積分を行列乗算に減らします。
xまず、方程式をから1までx、次に から0までで2 回積分できることは明らかです。したがって、次の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十分な大きさを選択しました。deltasayを選択した場合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または他の手段を使用して反復的に改善することもできます。この方法では、比較的高速な収束が得られ、膨大な計算を構築して解決する必要なく、必要な精度に到達できます。行列。