Solving Integral Equations - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: HP Prime Software Library (/forum-15.html) +--- Thread: Solving Integral Equations (/thread-14770.html) |
Solving Integral Equations - Eddie W. Shore - 04-03-2020 12:16 PM 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:
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/hp-prime-solving-integral-equations.html RE: Solving Integral Equations - Albert Chan - 04-03-2020 02:52 PM 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 RE: Solving Integral Equations - Albert Chan - 04-04-2020 05:58 PM 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} RE: Solving Integral Equations - Eddie W. Shore - 04-08-2020 02:32 PM I didn't realize I did this twice. RE: Solving Integral Equations - peacecalc - 11-03-2023 02:28 PM 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. |