Numerical methods, approximation in the sense of the least squares to find Polynomial

Question

Given these table of points x = [-1 -1/2 0 1/2 1] and f(x) = [2 1 0 1 2]. Find the approximation polynomial in the sense of the least square with weight function $w(x) = 1.$

My attempt:

So we take the canonical base $1, x, x^2$ and thus we have the polynomial:

$$P(x) = c_0 + c_1x+c_2x^2$$

Our goal is to find $c_0, c_1, c_2$.

We write the system: $$\left\{\begin{array}{l}c_0<1, 1>+c_1 <1, x> + c_2 <1, x^2>= <1, f>\\c_0<x,1>+c_1<x,x>+c_2<x,x^2>=<x,f>\\c_0<x^2,1>+c_1<x^2,x> +c_2<x^2,x^2>=<x^2,f>\end{array}\right.$$

Now we just calculate all the scalar products and find $c_0,c_1,c_2.$ My question is: how do I calculate the scalar product between 2 functions with discontinuous points? Do I just sum up the values of $x$ or $f$?And why is that? What is the intuition behind that?

Answer 1

Using the scalar product $$\langle f,g\rangle=\sum_{k=1}^5 w(x_k)f(x_k)g(x_k)$$ on the sample space gives the linear system $$ \pmatrix{5&0&2.5\\0&2.5&0\\2.5&0&2.125} \pmatrix{c_0\\c_1\\c_2} = \pmatrix{6\\0\\4.5} $$ leading to $c_1=0$. Eliminating $c_0$ gives $1.75c_2=3\implies c_2=\frac{12}7=1.7142857..$ and lastly $5c_0=6-\frac{30}7=\frac{12}7\implies c_0=\frac{12}{35}=0.342857..$.

Answer 2

I don't think you want to make this system too hard. There is a much more easy way to solve the least squares. I hope this answer is also applicable to your problem, since I'm not following your system. At least the final polynomial is in this answer, so that can be seen as a hint or a direction.

The least squares formula is given by minimizing \begin{equation} \sum_{i=1}^n y_i - f(x_i,\beta) \end{equation} which can be done by differentation and solving for 0, which yields, in matrix form, the matrix of coefficients $c$: \begin{equation} c = (A^TA)^{-1}A^Ty \end{equation} I assume you are familiar with this. Now, let's see how to solve this numerically. I programmed it in R, but I hope you can see how to write this program in any langauge.

library(plotly)
library(dplyr)

Non_Lin_Regressie <- function(x,y){
  # Solving the least squares method
  Pol_Matrix          <-    matrix(c(x^2,x,rep(1,length(x))),ncol=3)
  LHS_pol             <-    t(Pol_Matrix) %*% Pol_Matrix            # Computes the LHS
  RHS_pol             <-    t(Pol_Matrix) %*% y                     # Computes the RHS
  x_sol_pol           <-    solve(LHS_pol,RHS_pol)                  # Solves the equation
  return(x_sol_pol)
}

x <- c(-1,-1/2,0,1/2,1)
y<- c(2,1,0,1,2)
sol <- Non_Lin_Regressie(x,y)
func <- sol[1]*x^2 + sol[2]*x+sol[3]
plot_ly(x=~x, mode= 'lines') %>%
  add_trace(y = ~func, line = list(shape = "spline"), name= 'solved by LS') %>%
  add_trace(y=~y, name= 'normal plot')

Which yields the following coefficients matrix

          [,1]
[1,] 1.7142857
[2,] 0.0000000
[3,] 0.3428571

And hence the following polynomial \begin{equation} 1.71x^2+0.34 \end{equation} We can also plot both the given data and our fitted polynomial