[Python] Discrepancy in Python and Maple SVR Results

"I'm working on solving integral equations using the Least Square SVR method. I'm supposed to convert the Maple code provided in the link below into Python. I've written my Python code as shown below, but no matter what I do, I'm not getting the exact results. Could you please tell me how I can make my code produce the same results as Maple?

https://github.com/alirezaafzalaghaei/LSSVR-FIE/blob/master/paper-examples/example-4/CLS-SVR-dual.mw

import numpy as np
from scipy.integrate import quad
from scipy.special import legendre
from scipy.optimize import minimize
import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.stats import qmc
from scipy.integrate import quad
from scipy.integrate import dblquad
import time
from tensorflow.keras import regularizers
from tensorflow.keras.losses import MeanSquaredError

# Set precision
np.set_printoptions(precision=15)

# Define integration function
def Quad(f, a, b):
# Note: This is a simplified version, you may need to adjust it based on the
specific quadrature method you want to use.
return dblquad(f, a, b,a,b)[0]

# Define interval
a, b = 0, 1

# Define exact function
def exact(x,y):
return 1/((x+y+1)**2)

# Define function f
def f(x,y):
return 1/((x+y+1)**2)-(x/(6*(8+y)))

# Define function k
def k(x, y,t,s):
return (x/((8+y)*(1+t+s)))

# Example usage (replace with your specific use case)
result_quad = Quad(f, a, b)
print(result_quad)
import numpy as np
from scipy.special import legendre
from scipy.optimize import fsolve

def shift(x):
return (2 * x - a - b) / (b - a)

gamma = 10 ** 8
M = 4
N = M + 1

a, b = 0, 1 # Define a and b

# Calculate the roots of the Legendre polynomial
train1 = fsolve(lambda x: legendre(N)(shift(x)), np.linspace(0, 1, N))
train2 = fsolve(lambda y: legendre(N)(shift(y)), np.linspace(0, 1, N))

# Create the N x 2 matrix using NumPy's meshgrid and stacking
X, Y = np.meshgrid(train1, train2)
train = np.stack((X.flatten(), Y.flatten()), axis=1)

print(train)
print(train.shape) # Output the shape to confirm
phi = []
L_phi = []
for kindx in range((M + 1) * (M + 1)):
i = kindx // (M + 1) + 1
j = kindx % (M + 1) + 1
#print(i,j)
# Create the Legendre polynomial of degree i
p = legendre(i)
# Define phi using a nested function to create a new scope
def make_phi(p=p): # Capture p in the nested function's argument
return lambda x: p(shift(x))
phi.append(make_phi())
L_phi.append(lambda x,y: phi(x)*phi[j](y) - dblquad(lambda t,s: k(x,y, t,s)
* phi(t)*phi[j](s), a, b,a,b)[0])

train = np.array([train1, train2]).T # Assuming train1 and train2 are lists or
arrays
A = np.zeros((M+1, M+1))
for m in range(M+1):
for n in range(M+1):
A[m, n] = L_phi[m](train[n, 0], train[n, 1])
print(A)
Omega = np.dot(A.T, A) + np.identity(M+1) / gamma
GAMMA = np.array([f(train1,train2) for i in range(M+1)]).reshape(M+1, 1)
alpha = solve(Omega, GAMMA)
import numpy as np
l=[]
import numpy as np
from scipy.linalg import solve

# ... (other imports and functions) ...

def u_tilde(x, y):

# Evaluate only the first N Legendre polynomials in phi at x and y
# This aligns with the dimensions of A
P_x = np.array([phi(x) for i in range(M+1)]).reshape(-1, 1) # Changed here
P_y = np.array([phi(y) for i in range(M+1)]).reshape(-1, 1) # Changed here

# Calculate u˜(x) using element-wise multiplication and summation
u_tilde_x = alpha.T @ A.T @ (P_x * P_y) # Changed here

return u_tilde_x[0, 0] # Extract the scalar value from the result

# ... (rest of the code) ...# Extract the scalar value from the result


return result
# Example usage:
x_value = train1
y_value = train2
for k in range(M+1):
u_tilde_at_x = u_tilde(x_value[k],y_value[k])
l.append(u_tilde_at_x)
print(f"u˜({x_value[k]}) = {u_tilde_at_x}")
#u_tilde_at_x = u_tilde(x_value)
#print(f"u˜({x_value}) = {u_tilde_at_x}")

Continue reading...