Solving Integral Equations

[Eddie W. Shore]

The program INTEGRALSOLVE solves the following integral equation for x:

x
∫ f(X) dX - a = 0
0

using Newton's method.

Big X represents the variable of f(X) to be integrated while small x is the x that needs to be solved for.

Taking the derivative of the above integral using the Second Fundamental Theorem of Calculus:

d/dx [ ∫( f(X) dX from X=0 to X=x ) - a ]
= d/dx [ F(x) - F(0) - a ]
= d/dx [ F(x) ] - d/dx [ F(0) ] - d/dx [ a ]
= d/dx [ F(x) ]
= f(x)

F(X) is the anti-derivative of f(X). F(0) and a are numerical constants, hence the derivative of each evaluates to 0.

Newton's Method to solve for any function g(x) is:

x_n+1 = x_n - g(x_n) / g'(x_n)

Applying this to the equation, Newton's Method gives:

x_n+1 = x_n - [ ∫( f(X) dX from X=0 to X=x_n ) - a ] / f(x_n)

HP Prime Program INTEGRALSOLVE

Note: Enter f(X) as a string and use capital X. This program is designed to be use in Home mode.

EXPORT INTEGRALSOLVE(f,a,x)

Code:


BEGIN
// f(X) as a string, area, guess
// ∫(f(X) dX,0,x) = a
// EWS 2019-07-26
// uses Function app
LOCAL x1,x2,s,i,w;
F0:=f;
s:=0;
x1:=x;
WHILE s==0 DO
i:=AREA(F0,0,x1)-a;
w:=F0(x1);
x2:=x1-i/w;
IF ABS(x1-x2)<1ᴇ−12 THEN
s:=1;
ELSE
x1:=x2;
END;
END;

RETURN approx(x2);
END;

Examples

Radians angle mode is set.

Example 1:

Solve for x:

x
∫ sin(X) dX = 0.75
0

Initial guess: 1

Result: x ≈ 1.31811607165

Example 2:

Solve for x:

x
∫ e^(X^2) dX = 0.95
0

Initial guess: 2

Result: x ≈ 0.768032819934

Blog link: http://edspi31415.blogspot.com/2020/04/h...tions.html

[Albert Chan]

Numerically getting F(x) is expensive, required possibly many f(x) calls.
We can use third order iterative formulas instead.

Redo your 2nd example (Mathematica code):

In[1]:= f[x_] := Exp[x^2]
In[2]:= F[x_] := NIntegrate[f[t], {t, 0, x}] - 0.95

In[3]:= order2[x_] := x - F[x]/f[x]         (* newton's method *)
In[4]:= NestList[order2, 2., 8]
Out[4]= {2., 1.71606, 1.39774, 1.07527, 0.84433, 0.772609, 0.768049, 0.768033, 0.768033}

In[5]:= order3[x_] := order3[x, F[x], f[x]]
In[6]:= order3[x_, Fx_, fx_] := x - Fx/mean[fx, f[x-Fx/fx]]

In[7]:= mean[a_, b_] := (a + b)/2         (* algorithm 7 *)
In[8]:= NestList[order3, 2., 6]
Out[8]= {2., 1.57877, 1.1098, 0.803199, 0.768074, 0.768033, 0.768033}

In[9]:= mean[a_, b_] := 2/(1/a + 1/b)   (* algorithm 9 *)
In[10]:= NestList[order3, 2., 5]
Out[10]= {2., 1.45024, 0.908868, 0.769107, 0.768033, 0.768033}

In[11]:= Last[%] // CForm
Out[11]//CForm= 0.768032819933785

[Albert Chan]

I noticed you had post this same thread earlier, on Aug 24, 2019
Back then, I suggested reusing previous integral calculations.

Combined with previous post Algorithm 9, we have:

In[1]:= f[x_] := Exp[x^2]

In[2]:= order3[{x_, x0_, c_}] := Module[{t, Fx, fx},
               Fx = NIntegrate[f[t], {t, x0, x}] + c;
               fx = f[x];
               {x - Fx/mean[fx, f[x - Fx/fx]], x, Fx}
            ]

In[3]:= mean[a_, b_] := 2/(1/a + 1/b) (* algorithm 9 *)

In[4]:= First /@ NestList[order3, {2, 0, -0.95}, 5]

Out[4]= {2, 1.45024, 0.908868, 0.769107, 0.768033, 0.768033}

[Eddie W. Shore]

I didn't realize I did this twice.

[peacecalc]

Hello Eddie,

thank you very much for your little program.

I added a loop counter for stopping the program after 100 integrations.
For controlling I added some PRINT() commands (and I was astohnished about how fast the the precision is reached).

And I used a another technique for integrating (with cas command "int") and instead of always integrating from the lower (fixed) value to the new upper value, the integration is cutted in parts and added up:

\[ \displaystyle\int_{xo}^{xn} F1(X) dX = \int_{xo}^{x1} F1(X) dX + \int_{x1}^{x2} F1(X) dX
+ ... + \int_{xn-1}^{xn} F1(X) dX
\]

With the hope, if xn-1 and xn converge, the numerical part of integration becomes faster and faster.