Post Reply 
Solving Integral Equations
04-03-2020, 12:16 PM (This post was last modified: 04-03-2020 12:16 PM by Eddie W. Shore.)
Post: #1
Solving Integral Equations
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
Visit this user's website Find all posts by this user
Quote this message in a reply
04-03-2020, 02:52 PM (This post was last modified: 04-03-2020 03:10 PM by Albert Chan.)
Post: #2
RE: Solving Integral Equations
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
Find all posts by this user
Quote this message in a reply
04-04-2020, 05:58 PM
Post: #3
RE: Solving Integral Equations
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}
Find all posts by this user
Quote this message in a reply
04-08-2020, 02:32 PM
Post: #4
RE: Solving Integral Equations
I didn't realize I did this twice.
Visit this user's website Find all posts by this user
Quote this message in a reply
11-03-2023, 02:28 PM
Post: #5
RE: Solving Integral Equations
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.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)