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
Post Reply 




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