#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 19 15:40:37 2021
Least square minimization of error, using linear elements to reduce
the error of f(x)-fhat(x), where fhat is the sum of a_i * phi_i, i=0 to n
We minimize by finding the zero of the error with respect to the a_i
This version simplifies f(x) by making it merely f(x)=sin(x)
This is the 4th and final program in fitting a piece-wise linear function
to a known function.
@author: neil
"""

# compares the minimization of the R-squared error with the simpler method of reducing
# the integral of the linear error to zero
import numpy as np
import matplotlib.pyplot as plt

#xi = [0.5,1,1.5,2,3,4,4.75,5.5,6,7]     # (testing)list of element boundaries
xi = [0.5,1.5,4,5.5,7]                   # list of element boundaries
xpoints = np.linspace(xi[0],xi[-1],100)  # for plotting, make a list of points
f = np.sin(xpoints)    # this is the 'f' function we are approximating

fig = plt.figure()           # plot the 'f' function in red, leave the axis to plot 'fhat'
ax  = fig.add_subplot(111)  
ax.plot(xpoints,f,"r-")
ax.set_ylabel("f(x)-red, fhat(x)-black (minimize residual), fhat(x)-blue (error sum)")
ax.set_title("Finite Element curve fitting, least squares and linear error")
ax.set_xlabel("x-coordinate, showing element boundaries")

for i in range(len(xi)):              # do each element in sequence
    ax.plot([xi[i],xi[i]],[-1.5,1.5],"g:") # add the element edges        

# *********** now make minima of R**2, this is the least square error method
A = np.zeros([len(xi),len(xi)])  # set up A and C matric/vectors for linear algebra
C = np.zeros(len(xi))

for k in range(1,len(xi)-1):
    le0 = xi[k] - xi[k-1] # element lengths, le0 is the element in between a[k] and a[k-1]
    le1 = xi[k+1] - xi[k]
    A[k,k-1] = le0/6      # integral of phi_i-1 * phi_i
    A[k,k]   = le0/3 + le1/3
    A[k,k+1] = le1/6
    # the only difficult integral is integral of f(x) * phi_i
    C[k]     = (np.sin(xi[k])-np.sin(xi[k-1]))/le0 - (np.sin(xi[k+1])-np.sin(xi[k]))/le1

A[0,0]   = 1   # Make the a_i at the ends be exact
A[-1,-1] = 1
C[0]     = np.sin(xi[0])
C[-1]    = np.sin(xi[-1])

a_i = np.linalg.solve(A,C)
ax.plot(xi,a_i,'k-d')

# *************** now do the integral method of making the sum of the errors go to zero
# this is the integral of the discete sum method of version 2
A = np.zeros([len(xi),len(xi)])
C = np.zeros(len(xi))
C[0] = np.sin(xi[0])
A[0,0] = 1

for i in range(1,len(xi)):              # do each element in sequence
    x1 = xi[i-1]
    x2 = xi[i]
    # integral of 'f' from x1 to x2
    If = -np.cos(x2)+np.cos(x1)
    b = (x2-x1)/2
    A[i,i]   = b
    A[i,i-1] = b
    C[i] = If
    
a = np.linalg.solve(A,C)    
ax.plot(xi,a,"b-d")