c99e3b99-ac50-407c-9f66-e14bb2ed6d9d 34 4 5 false false 0 decimal

&[DATE] &[TIME] - &[FILENAME]

Integration using ode solvers

Integrals can be evaluated with an appropiate ode solver. For example:

x φ x x cos * :

a b 12 mat 0 π 12 mat :

x φ x a b int 2 -

x y D x φ 11 mat : M 011 mat a b 100 x y D rkfixed : M M rows 2 el 31 line 2.04951252 -

Notice that Rkadapt fails, but not this other from alglib

x y D x φ 11 mat : M 011 mat a b 100 x y D al_rkckadapt : M M rows 2 el 31 line 2 -

The ode solver

We can define our own integration procedure, cleaning the ode solver method for improve it's performance, This is a rude Runge Kuta 2-3 solver

u.0 t.1 t.2 N E RK23 u u.0 : t t.1 : h.o t.2 t.1 - N / eval : h h.o : 14 mat t k K h t D * eval h k * eval try 11 line : t t.2 < k.1 t u K : k.2 t h + k.1 K : k.3 t 0.5 h * + k.2 K : d k.1 2 k.3 * - k.2 + 3 / abs eval : (q E u abs 112 mat max * eval : (≤ t t h + : u u k.1 4 k.3 * + k.2 + 6 / + eval : 21 line if d 0 ≠ h h.o 0.9 h * q d / 3 nthroot * 12 mat min eval : if 51 line while u 41 line :

Thus

x D x φ : 0 a b 100106 -^RK23 21 line 2.00113492 -

Double Integrals

SMath can Integrate

x y f x 2^y 2^+ :

over the region defined between

x a e x^: x b 2 x * 2 + : 21 sys

x.1 x.2 21 mat x b x a - x solve sort : 0.768 - 1.6783 21 mat

x a x b x.1 x.2 2 1 mat x.1 x.2 2 1 mat a vectorize o augment 3 1 sys

t.0 0 time :

Numerical optimization x y f y x a x b int x x.1 x.2 int evalf maple 18.57686882

I.0 x y f y x a x b int x x.1 x.2 int eval : 18.57686864

0 time t.0 - 1.891 s *

Define some utilities

t.0 Err\0026\T Δt 0 time t.0 - : t.0 0 time : 12 mat I I.0 - abs Δt 21 mat 21 line :

ε.x 103 -^:

ε.y 103 -^:

n.x 15 :

n.y 15 :

Now, we try to extend the above method for solve this problem parametrizing the ode solver from the SMath plugins ...

2 f x.1 x.2 y.1 y.2 n.x n.y DInt a b n RK M 011 mat a b n x y D al_rkckadapt : M M rows 2 el eval 21 line : y u D x v D x y f 11 mat : x.1 x.2 n.x RK 21 line : y.1 y.2 n.y RK 31 line :

... or using our own ode solver ...

2 f x.1 x.2 y.1 y.2 n.x n.y ε.x ε.y DInt y D x D x y f : 0 x.1 x.2 n.x ε.x RK23 21 line : 0 y.1 y.2 n.y ε.y RK23 21 line :

... and, as first approach, using integration over a rectangle bigger than the true region of integration and making zero f outside the region of integration

x y g x y f y x b ≤ (* y x a ≥ (* :

y.1 y.2 21 mat 0621 mat :

t.0 0 time :

I x y g y y.1 y.2 int x x.1 x.2 int : 18.5832591

t.0 Err\0026\T 0.0064 5.35 s * 21 mat

I x y g x.1 x.2 y.1 y.2 n.x n.y DInt : #

Convergent only for rkafixed ode solver.

I x y g x.1 x.2 y.1 y.2 n.x n.y ε.x ε.y DInt : 18.45199634

t.0 Err\0026\T 0.1249 71.53 s * 21 mat

As better approach, using a change of variables

y y.1 u y.2 y.1 - (* + ≡

dy y.2 y.1 - (du * ≡

we can integrate from 0 to 1 for y, and isn't necessary to guess the limits of integration, as in the first method.

x u h x x a u x b x a - (* + f x b x a - (* :

I x y h y 01 int x x.1 x.2 int : 18.56839179

t.0 Err\0026\T 0.0085 10.453 s * 21 mat

I x y h x.1 x.2 01 n.x n.y DInt : 18.57686894

t.0 Err\0026\T 3.0316 107 -^* 7.387 s * 21 mat

I x y h x.1 x.2 01 n.x n.y ε.x ε.y DInt : 18.4681218

t.0 Err\0026\T 0.1087 2.139 s * 21 mat

Alvaro