What's wrong with this code for solving integral equations in python

Question

I have the following integral equation:

I'm trying to solve the integral equation to see whether f(u) is predicted as cos(2u). Since we know that the solution to the integral equation is cos(2u), we can approximate the integral from 0 to infinity to the limits 0 to say, 5 if we make the value of the integral from 5 to infinity negligible, and this can be done by choosing t to be small. I have chosen 100 evaluation points for the integral between 0 to 5, and this implies that i am solving for 100 values of f(u). Since i need to solve for 100 values of f(u), I need to generate 100 equations, and thus need 100 values of time t. I choose 100 values for time t between 1 and 1.3 since this will ensure that the integral is negligible for values of 5 and beyond. The following is the scipy code for doing this:

from scipy import*
from matplotlib.pyplot import*


Nt_samples=100  #100 evaluation points for the time t
t=linspace(1.0,1.3,100)
number_eval_points=100  #100 evaluation points for u 

eval_points=linspace(0.005,5,number_eval_points)
delta=eval_points[1]-eval_points[0]
R=zeros(100,1)
R=0.5*sqrt(2*3.14)*t*exp(-2*t*t)
A=zeros((Nt_samples,number_eval_points))    

for i in range(100):
    for j in range(100):
        A[i,j]=delta*exp(-(eval_points[j]*eval_points[j])/(2*(t[i]*t[i])))


Z=cos(2*eval_points)
Fu=dot(linalg.inv(A),R)
plot(eval_points,Fu,eval_points,Z)

Somehow, my results for f(u) are a far cry from cos(2u). In fact, they look like a lot of random noise and follow no pattern at all! Also, the magnitude of f(u) is very large. I have tried playing around with the number of evaluation points and the value of t, but I have no luck.

Is there anything wrong with the code/setting of parameters/logic?

Thanks a million!

Answer 1

This is not a programming question, but here goes: integral equations are quite often numerically badly conditioned.

Indeed, in your case,

u, s, vh = linalg.svd(A)
print(s.max()/s.min())
# -> 4.03758642411e+16

This is the condition number, and it's huge. The matrix A is nearly singular, so there will be a large error in the solution.

Googling for "Tikhonov regularization" should get you started on how people work around issues like this. Solving integral equations is a mature field in mathematics, so googling should help you a lot here.

A quick regularization is replacing linalg.inv(A) with

linalg.pinv(A, 1e-8)

This gives something more cosine looking. The magic value 1e-8 depends on the integral kernel, but when things are about rounding error, good values to try are around sqrt(finfo(float).eps) which is saying that you trust half of the ~15 digits that floating point numbers have.

Also, a change of variables may be a better idea than truncating the range of u to [0, 5]. The numerical issue in your case is possibly related to the exponentially decreasing weight factor, so getting rid of that via a change of variables may also make the condition number better.