我可以建议一种将您的方程简化为积分方程的方法，该积分方程可以通过用矩阵近似其内核来进行数值求解，从而将积分减少为矩阵乘法。
首先，很明显，方程可以在 上积分两次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足够大，以便错误可见。如果您选择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或其他方法迭代改进——这样你将有一个相对较快的收敛，将能够达到所需的精度，而不需要构造和求解巨大的矩阵。