cms.interpolate.interpolate_lagrange のソースコード

import sys
import csv
import numpy as np
from math import exp, sqrt, sin, cos, pi
import matplotlib.pyplot as plt


#===================
# constants
#===================
xmin = -2.0
xmax =  2.0
norder = 4
xomin =  0.0
xomax =  2.2
noutput = 101


if __name__ == "__main__":
    if len(sys.argv) >= 2:
        norder = int(sys.argv[1])

outfile  = "interpolate_lagrange.csv"


def func(x):
    return 1.0 / (1.0 + x*x)




def interpolate_lagrange(x, y, x0):
    ndata = len(x)
    if ndata < 3:
        print("Error: ndata [{}] must be larger than 2".format(ndata))
        exit()
    w = 0.0;
    for i in range(0, ndata):
        w1 = 1.0
        for j in range(0, ndata):
            if i != j:
                w1 *= (x0 - x[j]) / (x[i] - x[j])
        w += w1 * y[i];
    return w


#===================
# main routine
#===================


def main():
    global xmin, xmax, n

    x = []
    y = []
    xstep = (xmax - xmin) / norder
    for i in range(0, norder+1):
        x.append(xmin + i * xstep)
        y.append(func(x[i]))

    xf = []
    yc = []
    yf = []
    xostep = (xomax - xomin) / (noutput - 1)
    for i in range(0, noutput):
        xf.append(xomin + i * xostep)
        yc.append(func(xf[i]))
        yf.append(interpolate_lagrange(x, y, xf[i]))

    print("{:^14}  {:^14}  {:^14}".format('x', 'f', 'f(fit)'), end = '')

    print("Write to [{}]".format(outfile))
    f = open(outfile, 'w')
    fout = csv.writer(f, lineterminator='\n')
    fout.writerow(['x', 'f', 'f(fit)'])
    for i in range(0, noutput):
        print("{:>14.4f}  {:>14.4f}  {:>14.4f}".format(xf[i], yc[i], yf[i]))
        fout.writerow([xf[i], yc[i], yf[i]])
    f.close()

    print("")

# create figure object and axes list
    fig, ax = plt.subplots(1, 1)
    plt.title("Runge's phenomenon: order {}".format(norder))
    plt.xlabel('x')
    plt.ylabel('1 / (1 + x^2)')
    plt.plot(xf, yc, label = 'exact')
    plt.plot(xf, yf, label = 'fit')
    plt.legend()
    plt.show()



if __name__ == '__main__':
    main()