Oyama zemi ex1 linear interpolation

In [1]:

    

;cat linear_interp.jl


    

    




function li_ikegami(grid, vals)
    n = length(grid)
    m = length(vals)

    if n != m
        print("grid size and value siza is not equal")
    end

    mat = [grid vals]
    mat = sortrows(mat, by = x->(x[1]))
    mat = transpose(mat)

    grid = mat[1, :]
    vals = mat[2, :]

    slope = Array(Float64, n-1)
    for i in 1: (n-1)
        temp_slope = (vals[i] - vals[i+1]) / (grid[i] - grid[i+1])
        slope[i] = temp_slope
    end
    function approx(x_val)

        if x_val > maximum(grid)
            return print("too big value ...")
        end

        if x_val < minimum(grid)
            return print("too small value ...")
        end

        interval = 1
        for j in 1:(n-1)
            if x_val > grid[j+1]
                interval += 1
            else break
            end
        end
        return slope[interval] * (x_val - grid[interval]) + vals[interval]
    end

    function approx{T<:Real}(x::AbstractVector{T})
        n = length(x)
        out = Array(Float64, n)
        for i in 1:n
            out[i] = approx(x[i])
        end
        return out
    end

    return approx
end

In [2]:

    

include("linear_interp.jl")


    

    
Out[2]:





li_ikegami (generic function with 1 method)

In [3]:

    

grid = [1, 2]
vals = [2, 0]
f = li_ikegami(grid, vals)

f(1.25)


    

    
Out[3]:





1.5

In [4]:

    

using Gadfly
set_default_plot_size(18cm, 12cm)


    

In [5]:

    

mesh = 200
grid = linspace(1,2,mesh)
vals = Array(Float64, mesh)
for i in 1: mesh
    vals[i] = f(grid[i])
end
a = 1
b = 2
plot([f], a, b)


    

    
Out[5]:

In [6]:

    

a = 0
b = 20
mesh = 15
s = Array(Float64, mesh)

grid = linspace(a,b,mesh)
for i in 1:mesh
    s[i] = sin(grid[i])
end

g = li_ikegami(grid, s)
plot([g, sin], a,b)


    

    
Out[6]:

In [7]:

    

a = 0
b = 20
mesh = 15
s = Array(Float64, mesh)

grid = linspace(a,b,mesh)
for i in 1:mesh
    s[i] = sin(grid[i])
end

g = li_ikegami(grid, s)
function residual(x)
    return sin(x) - g(x)
end
plot([residual], a,b)


    

    
Out[7]:

In [9]:

    

function sum_resi(mesh)   
    a = 0
    b = 20
    mesh_2 = 200
    
    s = Array(Float64, mesh)
    grid = linspace(a,b,mesh)
    
    for i in 1:mesh
        s[i] = sin(grid[i])
    end
    g = li_ikegami(grid, s)
    
    function residual_sq(x)
        return (sin(x) - g(x))^2
    end
    
    grid_2 = linspace(a,b,mesh_2)
    w = Array(Float64, mesh_2)
    
    for i in 1:mesh_2
        w[i] = residual_sq(grid_2[i])
    end 
    return sum(w)
end

max = 20
x_line = collect(2:max)
q = Array(Float64, max-1)

for i in x_line
    q[i-1] = sum_resi(i)
end
plot(x = x_line, y = q, Guide.xlabel("number of mesh"), Guide.ylabel("squared residual"), Guide.title("Decreasing Squared Residual"))


    

    
Out[9]:

In [ ]: