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
 Eddie W. Shore Senior Member Posts: 1,122 Joined: Dec 2013
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

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

04-03-2020, 02:52 PM (This post was last modified: 04-03-2020 03:10 PM by Albert Chan.)
Post: #2
 Albert Chan Senior Member Posts: 970 Joined: Jul 2018
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
04-04-2020, 05:58 PM
Post: #3
 Albert Chan Senior Member Posts: 970 Joined: Jul 2018
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}
04-08-2020, 02:32 PM
Post: #4
 Eddie W. Shore Senior Member Posts: 1,122 Joined: Dec 2013
RE: Solving Integral Equations
I didn't realize I did this twice.
 « Next Oldest | Next Newest »

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