n=50
a=-5;b=5;x=0*1:n

for ( k in 1:n){  # Nós de Chebyshev reescalados.
x[k]=(a+b)/2+(b-a)*cos((2*k-1)/(2*n)*pi )/2}

y=1/(1+x^2)

n=length(x)
plot(x,y,col="red")

x=x/5      # Mudança de varivál paa melhorar condicionamento do sistema
plot(x,y,col="red")

# Agora cálculos mostram que a função g(u)= |Au-y|^2 / 2
# em que A é a matriz e é um vetor, como montamos no que segue.
n=length(x)
A=matrix(0,n,n); 
for ( j in 1:n){A[,j]=x^(j-1)}


u=0*x; v=y
B=t(A)%*%A; b=t(A)%*%v           # Adaptando o sistema para A*Au=A*v, ou Bu=b, em que B=A*A é positiva definida e b=A*v.
                                 # Início do método
w=B%*%u-b                        # w0
d=w                              # d1
h= t(w)%*%w /(t(B%*%w) %*%w)     # Escolha de h0 ou alpha0
u=u-h[1]*w                       # u1
w=B%*%u-b                        # w1

for ( j in 1:100*n){                
bt= t(B%*%d) %*%w/(t(B%*%d) %*%d) # Escolha de beta
d=w-bt[1]*d  
at= t(d) %*%w/(t(B%*%d) %*%d)    # Escolha de alpha
u=u-at[1]*d                   
w=B%*%u-b  
}
u                                # Aproximação para u(t), sendo t a soma dos passos h.
Erro1=A%*%u-v 
print("Residuo 1")
Erro1   

poli<-function(s){
    p=u[1]
    for ( i in 2:n){p=p+u[i]*(s/5)^(i-1)}
    return(p)
} 

curve(poli,-5,5,col="blue",ylim=c(0,1))
points(5*x,y,col="red")        # Teste de ajuste.


Erro<-function(s){1/(1+s^2)-poli(s)}
curve(Erro,-5,5)

# Se quiser pode aumentar t e n ou escolher outra condição inicial 
# para se convencer do resultado obtido.

#----------------------------------------------------------------------------
# Agora cálculos mostram que a função g(u)= |f(u)|_1, com f(u)=Au-y
# 

f<-function(u){A%*%u-y}

ll=10^(-15)

sf<-function(i){
    p=0*A[1,]
    q=A[i,]%*%x-y[i]
    if (q[1]>ll){p=A[i,]}
    if (-q[1]>ll){p=-A[i,]}
    return(p)
}
SGradf<-function(u){
    p=0*A[1,]
    for (i in 1:n){p=p+sf(i)}
    return(p)
}

print("Teste subgrad")
SGradf(y)

# Método dos subgradientea
vv=u
for ( j in 1:(10000)){ 
w=SGradf(u) 
mw=t(w)%*%w
u=u-w/sqrt(mw[1]) 
if (sum(abs(f(u)))<sum(abs(f(vv)))){vv=u}
}
vv                                # Aproximação para u(t), sendo t a soma dos passos h.
Erro2=A%*%vv-v 
print("Residuo 2")
Erro2    
u=vv

polis<-function(s){
    p=u[1]
    for ( i in 2:n){p=p+u[i]*(s/5)^(i-1)}
    return(p)
} 

curve(polis,-5,5,col="blue",ylim=c(0,1))
points(5*x,y,col="red")        # Teste de ajuste.


Erro<-function(s){1/(1+s^2)-poli(s)}
curve(Erro,-5,5)

# Se quiser pode aumentar t e n ou escolher outra condição inicial 
# para se convencer do resultado obtido.

abs(Erro2-Erro1)