Hello everyone!

I am Ludovico. I have been accepted to the Julia season of contribution, and during the community bonding time I started to tackle the first part of my proposal. I am building a new package called “Surrogates.jl”. My goal is to implement methods to approximate solution of Differential equations.

I started developing two methods: “Kriging” and “Response surfaces”. They both deal with radial basis functions and polynomials. As of right now, I just have the code for the one dimensional case, but the general case is not much more difficult. I am pretty sure that the code I have written is not “Julia-like” at all, so I would love to hear all the suggestions you have for it, if you do not mind.

I guess I should just post the first function, because I think that all the suggestions would work for the second one as well.

Below you can take a look at the code:

```
function Radial_1D(x,y,a,b,kind::String,lambda = 0)
#=
(x,y) set of nodes
(a,b) interval
kind is type of radial basis function
lambda is optional parameter with kind == multiquadric
=#
#1D Chebyshev is suggested
Chebyshev(x,k) = cos(k*acos(-1 + 2/(b-a)*(x-a)))
#Type of Radial basis function
#q is the number of polynomials in the basis
#The numbers are suggested by papers
if lambda === nothing
if kind == "linear"
q = 1
phi = z -> abs(z)
elseif kind == "cubic"
q = 2
phi = z -> abs(z)^3
elseif kind == "thinplate"
q = 2
phi = z -> z^2*log(abs(z))
else
error("Wrong type")
end
else
if kind == "multiquadric"
q = 1
phi = z -> sqrt(abs(z)^2 + lambda^2)
else
error("Wrong type")
end
end
if length(x) != length(y)
error("Data length does not match")
end
n = length(x)
#Find coefficients for both radial basis functions and polynomial terms
size = n+q
D = zeros(Float32, size, size)
d = zeros(Float32,size)
#In this array I have in the first n entries the coefficient of the radial
#basis function, in the last q term the coefficients for polynomial terms
coeff = zeros(Float32,size)
#Matrix made by 4 blocks:
#=
A | B
B^t | 0
A nxn, B nxq A,B symmetric so the matrix D is symmetric as well.
=#
for i = 1:n
d[i] = y[i]
for j = 1:n
D[i,j] = phi(x[i] - x[j])
end
for k = n+1:size
#Symmetric
D[i,k] = Chebyshev(x[i],k)
D[k,i] = D[i,k]
end
end
#Vector of size n + q containing in the first n terms the coefficients of
# the radial basis function and in the last q term the coefficient for the
# polynomials
return D\d
end
```