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,559 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: 2,498 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: 2,498 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,559 Joined: Dec 2013
RE: Solving Integral Equations
I didn't realize I did this twice.
11-03-2023, 02:28 PM
Post: #5
 peacecalc Member Posts: 202 Joined: Dec 2013
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.
 « Next Oldest | Next Newest »

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