<a href="https://colab.research.google.com/github/sergiorolnic/github-slideshow/blob/master/Untitled11.ipynb" target="_parent"></a>

# Esempi

## prime due es


```python
#Scrivere un programma Python per calcolare la distanza tra due punti

import math
p1=[4,0]
p2=[6,6]

distanza= math.sqrt((p1[0]-p2[0])**2+(p1[1]-p2[1])**2)


Creare una lista con 100 elementi uguali a 0.


lista=[0]*100
print(lista)

'''
Esercizio 20
Data una lista annidata a=[[2,1,3],[4,5,6],[7,8,9]], scrivere un programma Python che la trasformi 
nella lista semplice ottenuta con i suoi valori, cioè a1=[2,1,3,4,5,6,7,8,9]
'''


a=[[2,1,3],[4,5,6],[7,8,9]]

a1=a[0]+a[1]+a[2]
print(a1)

z = np.arange(0,30.5,0.5)
print("Numpy Array z ",z)
t = np.linspace(1,2,100)

ps= np.dot(a,b)
print("Prodotto scalare tra a e b ",ps)

D=np.diag(a)


a=1
c=1

k=np.arange(1,9)
b=10.0**k;
Delta= b**2-4*a*c
'''
La formula del calcolo della soluzione x1 può dare problemi numerici in base ai valori di b:
quando b assume un valore molto elevato b**2-4ac è approssimabile a b**2, quindi nel calcolo
di x1 si può veriificare il fenomento di cancellazione dovuto al calcolo -b+ sqrt(Delta)
in quanto sqrt(Delta) in questo caso è approssimabile a b.
Per la formula del calcolo della soluzione x2 questo problema non si verifica.
'''
x1=(-b+np.sqrt(Delta))/(2*a)
x2=(-b-np.sqrt(Delta))/(2*a)

'''
Usiamo una formula alternativa, nota in algebra, 
per il calcolo della soluzione x1 a partire dalla soluzione x2
'''

x1new=c/(a*x2)

------

```



# Malcondizionamento sistema al variare di x




```python
import numpy as np
import numpy.linalg as npl
import scipy.linalg  as spl
import matplotlib.pyplot as plt

A=np.array([[3,5],[3.01,5.01]])
b=np.array([10, 1])

x= spl.solve(A,b)  #soluzione del sistema: x=[-2255; 1355]

# perturbo il coeff. di x della seconda equazione:
dA=np.array([[0, 0], [0.01, 0]])


err_dati= npl.norm(dA,np.inf)/npl.norm(A,np.inf)  #errore percentuale=0.12%
print("Errore relativo sui dati Sistema 1 pert. in percentuale", err_dati*100,"%")
# soluzione del sistema perturbato:
x1=spl.solve(A+dA,b) 


err_rel_sol=npl.norm(x-x1,np.inf)/npl.norm(x,np.inf) #%errore percentuale=71.43%
print("Errore relativo sulla soluzione Sistema 1 pert. in percentuale", err_rel_sol*100,"%")
# quindi ad una piccola perturbazione sui dati corrisponde una grossa
# perturbazione sui risultati --> problema mal condizionato


"""
Esercizio 4

"""
import numpy as np
import numpy.linalg as npl
import scipy.linalg as spl

#Esempi di sistemi malcondizionati


A=np.array([[6, 63, 662.2],[63, 662.2, 6967.8],[662.2, 6967.8, 73393.5664]])
b=np.array([1.1, 2.33, 1.7])

KA= npl.cond(A,np.inf)

x=spl.solve(A,b)

#perturbare la matrice
 
A1=A.copy()
A1[0,0]=A[0,0]+0.01
x_per=spl.solve(A1,b)


#Errore relativo sui dati
err_dati=npl.norm(A-A1,np.inf)/npl.norm(A,np.inf)
print("Errore relativo sui dati  in percentuale ", err_dati*100,"%")

err_rel_sol=npl.norm(x_per-x,np.inf)/npl.norm(x,np.inf)
print("Errore relativo sulla soluzione  in percentuale ", err_rel_sol*100,"%")




```

    Errore relativo sui dati Sistema 1 pert. in percentuale 0.12468827930174566 %
    Errore relativo sulla soluzione Sistema 1 pert. in percentuale 71.42857142857439 %




# Zeri

## esami


```python
"""
Ese1
"""

import numpy as np
import matplotlib.pyplot as plt
from funzioni_zeri import newton_m, stima_ordine
import sympy as sym
from sympy.utilities.lambdify import lambdify
from scipy.optimize import fsolve

x=sym.symbols('x')
fx= x-1/3*sym.sqrt(30*x-25)
dfx=sym.diff(fx,x,1)
print(dfx)

x0=4
f=lambdify(x,fx,np)
alfa=fsolve(f,x0)
print("La funzione ha uno zero in ", alfa)
fp=lambdify(x,dfx,np)
fp_alfa=fp(alfa)
print("La derivata prima  in ", alfa, 'vale ',fp_alfa)
x=np.linspace(5/6,25/6,100)
plt.plot(x,f(x))
plt.plot(x,[0]*100)
plt.plot(alfa,0,'ro')
plt.show()


#In alfa=1.66666667 si annula sia la funzione che la sua derivata prima,
#la funzione ha in x=1.66666667 uno xero con molteplicita=2 
m=2
tolx=1e-12
tolf=1e-12
nmax=100

#Metodo iterativo che converge quadraticamente al alfa: metodo di Newton Modificato con m=2
x1,it,xk= newton_m(f,fp,x0,m,tolx,tolf,nmax)
#Verifico l'ordine di convergenza
ordine=stima_ordine(xk,it)

plt.plot(range(it),np.abs(xk))
plt.show()


#Il metodo non converge se scelgo come iterato iniziale x0=5/6,
#perchè la derivata prima in  5/6 diverge  va a -infinito
```


```python
# -*- coding: utf-8 -*-
"""
Created on Mon May 24 19:09:59 2021

@author: damia
"""

import numpy as np
import matplotlib.pyplot as plt
from funzioni_zeri import iterazione
import sympy as sym
from sympy.utilities.lambdify import lambdify
n=35
u1=np.zeros((n,),dtype=float)
u2=np.zeros((n,),dtype=float)
u3=np.zeros((n,),dtype=float)


for i in range(1,n+1):
    u1[i-1]=15*((3/5)**(i)+1)/(5*(3/5)**(i)+3)

u2[0]=4
for i in range(1,n):
    u2[i]=8-15/u2[i-1]


u3[0]=4
u3[1]=17/4
for i in range(3,n+1):
    u3[i-1]=108-815/u3[i-2]+1500/(u3[i-2]*u3[i-3])

plt.plot(range(n),u2)
plt.title('Formula 2')
plt.show()

plt.plot(range(n),u3)
plt.title('Formula 3')
plt.show()

err_rel2=np.abs(u2-u1)/np.abs(u1)

err_rel3=np.abs(u3-u1)/np.abs(u1)

plt.semilogy(range(n),err_rel2,range(n),err_rel3)
plt.legend(['Errore relativo formula 2', 'Errore relativo formula 3'])
plt.show()

g1=lambda x:  8-15/x

g2=lambda x:  108-815/x+ 1500/(x**2);

x=sym.symbols('x')

#Definisco funzione g1
g1x= 8-15/x
dg1x=sym.diff(g1x,x,1)
dg1=lambdify(x,dg1x,np)
g1=lambdify(x,g1x,np)


#Definisco funzione g2
g2x=108-815/x+1500/(x**2)
dg2x=sym.diff(g2x,x,1)
dg2=lambdify(x,dg2x,np)
g2=lambdify(x,g2x,np)

 
x=np.linspace(4,100 ,100)
plt.plot(x,g1(x))
plt.plot(x,x)
plt.legend(['g1(x)','y=x'])
plt.show()

#La g2 interseca la bisettrice in 2 punti, ha due punti fissi (x=5 ed 100) , ma la derivata prima della g2
# non soddisfa le ipotesi del teorema di convergenza locale in un intorno del primo punto fisso,5
plt.plot(x,g2(x))
plt.plot(x,x)
plt.legend(['g2(x)','y=x'])
plt.show()
 
tolx=1e-5
nmax=100
x0=4
x1,xk,it=iterazione(g1,x0,tolx,nmax)
x2,xk2,it2=iterazione(g2,x0,tolx,nmax)
print("Punto fisso della funzione g1 -->  ",x1)
print("Punto fisso della funzione g1 -->  ",x2)

#Visualizziamo la derivata prima di g1 in un intorno di 5, sono soddisfatte le iptesi del teorema di convergenza
#locale
xx=np.linspace(2,6,100)
plt.semilogy(xx,dg1(xx))
plt.plot([2,6],[1,1])
plt.plot([2,6],[-1,-1])
plt.legend(['derivata prima di g1 in un intorno di 5 ', 'y=1','y=-1'])
plt.show()
xx=np.linspace(2,6,100)
plt.semilogy(xx,dg2(xx))
plt.plot([2,6],[1,1])
plt.plot([2,6],[-1,-1])
plt.legend(['derivata prima di g2 in un intorno di 5 ', 'y=1','y=-1'])
 
    
plt.show()

xx=np.linspace(95,105,100)
plt.plot(xx,dg2(xx))
plt.plot([95,105],[1,1])
plt.plot([95,105],[-1,-1])
 
plt.legend(['Derivata prima  di g2 in un intorno di 100','y=1','y=-1'])
  


```


```python
import numpy as np
import sympy as sym
import funzioni_zeri
import matplotlib.pyplot as plt
from sympy.utilities.lambdify import lambdify
from scipy.optimize import fsolve
tolx=1e-7
nmax=1000


f= lambda x: np.tan(3/2*x)-2*np.cos(x)-x*(7-x)

#Utilizzo il metodo fsolve di scipy.optimize per calcolare lo zero alfa della funzione f,
#prende in input l'iterato iniziale x0

x0=0.0
alfa=fsolve(f, x0)
print("Lo zero della funzione e' ",alfa)
#Disegno: l'asse x e la funzione f valutata in un intervallo opportuno
xx=np.linspace(-1.0,1.0,100)
plt.plot(xx,0*xx,xx,f(xx),alfa,0,'ro')
plt.legend(['Funzione f', 'zero'])
plt.show()

'''
Definisco la funzione g in formato simbolico perchè poi utilizzo la funzione diff
di sympy per calcolare l'espressione analitica della derivata prima'
'''
x=sym.symbols('x')

#Considero la funzione g indicata dalla traccia del compito
gx=sym.tan(3/2*x)-2*sym.cos(x)-x*(6-x)

#Disegno la funzione g(x) e la bisettrice y=x
g=lambdify(x,gx,np)
plt.plot(xx,xx,xx,g(xx))
plt.title('funzione g(x) e y=x')
plt.show()

#Calcolo la derivata prima di gx espressione simbolica tramite la funzione diff del modulo sym 
dgx=sym.diff(gx,x,1)


dg=lambdify(x,dgx,np)

#Disegno la funzione dg(x) 

#Posso giustifcare la convergenza del procedimento iterativo guardando la derivata prima di g(x)
#in un intorno della soluzione: il metodo genera una successione di iterati convergenti alla radice alfa
# ed appartenenti a questo intorno se |g'(x)|< 1 in un intorno della soluzione

 
plt.plot(xx ,dg(xx ))
plt.plot(alfa,0,'ro')
#Disegno la retta y=1
plt.plot([-1,1],[1,1],'--')
#Disegno la retta y=-1
plt.plot([-1,1],[-1,-1],'--')
plt.title('funzione dg(x) proposta dalla traccia - Ipotesi per la convergenza non soddisfatte')
 
plt.legend(['Grafico derivata prima di g1 (x)', 'Zero', 'Retta y=1', 'Retta y=-1'])
plt.show()

#Dal grafico vedo che per la funzione g proposta dalla traccia
#non sono soddisfatte le ipotesi del teorema di convergenza locale


#Ricavo la funzione gx per la quale ci sia convergenza
gx1= (sym.tan(3.0/2.0*x)-2*sym.cos(x)+x**2)/7
#Disegno la funzione g(x) e la bisettrice y=x
g1=lambdify(x,gx1,np)
plt.plot(xx,xx,xx,g1(xx))
plt.title('funzione g1(x) ricavata e y=x')
plt.show()

#Calcolo la derivata prima di gx espressione simbolica tramite la funzione diff del modulo sym 
dgx1=sym.diff(gx1,x,1)


dg1=lambdify(x,dgx1,np)

#Disegno la funzione dg1(x) 

#Posso giustifcare la convergenza del procedimento iterativo guardando la derivata prima di g(x)
#in un intorno della soluzione: il metodo genera una successione di iterati convergenti alla radice alfa
# ed appartenenti a questo intorno se |g'(x)|< 1 in un intorno della soluzione

 
plt.plot(xx ,dg1(xx ))
plt.plot(alfa,0,'ro')
#Disegno la retta y=1
plt.plot([-1,1],[1,1],'--')
#Disegno la retta y=-1
plt.plot([-1,1],[-1,-1],'--')
plt.title('funzione dg1(x) Ricavata - Ipotesi di convergenza soddisfatte')
plt.legend(['Grafico derivata prima di g1 (x)', 'Zero', 'Retta y=1', 'Retta y=-1'])
plt.show()

'''
Dal grafico vedo che per la funzione g  ceh ho ricavato
soddisfatte le ipotesi del teorema di convergenza locale: esiste un intorno della soluzione per cui
|g'(x)|<1 

'''

#Utilizzo il metodo di iterazione funzionale per calcolare il punto fisso di g1
x1,it,xk=funzioni_zeri.iterazione(g1,x0,tolx,nmax)
print('iterazioni= {:d}, soluzione={:e} \n\n'.format(it,x1))

#Calcolo l'ordine del metodo
ordine_iter= funzioni_zeri.stima_ordine(xk,it)
#Essendo il metodo con ordine di convergenza lineare, la costante asintotica di convergenza è data
#da |g'(alfa)| dove alfa è la radice.

 
print("Iterazione it={:d}, ordine di convergenza {:e}".format(it,ordine_iter ))

plt.plot(range(it+1),xk)
 
```

## Newton, Newton modificato(molteplicità) e stima ordine


```python
import math
import numpy as np
import sympy as sym
import funzioni_zeri
import matplotlib.pyplot as plt
from sympy.utilities.lambdify import lambdify


tolx=1e-12;
tolf=1e-12;


x=sym.symbols('x')

fx=fx=x**3+x**2-33*x+63
dfx=sym.diff(fx,x,1)

#Trasformo in numeriche la funzione e la sua derivata
f=lambdify(x,fx,np)
df=lambdify(x,dfx,np)

#Disegno nello stesso grafico l'asse x e la funzione f valutata in un intervallo opportuno, [-10,10]
z=np.linspace(-10,10,100)
plt.plot(z,0*z,z,f(z),'r-')

nmax=500
x0=1
xNew,itNew,xkNew=funzioni_zeri.newton(f,df,x0,tolx,tolf,nmax)
print('X0= {:e} , zero Newton= {:e} con {:d} iterazioni \n'.format(x0,xNew,itNew))
ordine_New =funzioni_zeri.stima_ordine(xkNew,itNew)
print("Newton it={:d}, ordine di convergenza {:e}".format(itNew,ordine_New))

#Utilizzando il metodo di Newton modifica e ponendo m uguale alla molteplicità della radice
# si ottiene un metodo con ordine di convergenza 2

m=2
xNew_m,itNew_m,xkNew_m=funzioni_zeri.newton_m(f,df,x0,m,tolx,tolf,nmax)
print('X0= {:e} , zero Newton Mod= {:e} con {:d} iterazioni \n'.format(x0,xNew_m,itNew_m))
ordine_New_m =funzioni_zeri.stima_ordine(xkNew_m,itNew_m)
print("Newton Mod it={:d}, ordine di convergenza {:e}".format(itNew_m,ordine_New_m))
```

## Costate asintotica


```python
"""
esercizio 6
"""

import numpy as np
import sympy as sym
import funzioni_zeri
import matplotlib.pyplot as plt
from sympy.utilities.lambdify import lambdify

tolx=1.e-7
nmax=1000


f= lambda x: x**3+4*x**2-10

#Disegno: l'asse x e la funzione f valutata in un intervallo opportuno
xx=np.linspace(0.0,1.6,100)
plt.plot(xx,0*xx,xx,f(xx))
plt.show()


x0=1.5; 

x=sym.symbols('x')
#Definisco le possibili espressioni della funzione f(x)=0 nella forma x=g(x)
#Per come ricavare le diverse g(X) , analizzare file pdf su virtuale con spiegazione teorica dell'esercizio 6

#................................
#gx=sym.sqrt(10/(x+4));       #p=1, 0.127229401770925
gx=1/2*sym.sqrt(10-x**3)      #p=1, C=0.511961226874885
#gx=(10+x)/(x**2+4*x+1);         #p=1, C=0.983645643784931
#gx=sym.sqrt(10/x-4*x)         # non converge

#Disegno la funzione g(x) e la bisettrice y=x
g=lambdify(x,gx,np)
plt.plot(xx,xx,xx,g(xx))
plt.title('funzione g(x) e y=x')
plt.show()

#Calcolo la derivata prima di gx espressione simbolica tramite la funzione diff del modulo sym 
dgx=sym.diff(gx,x,1)


dg=lambdify(x,dgx,np)

#Disegno la funzione dg(x) 

#Posso giustifcare la convergenza del procedimento iterativo guardando la derivata prima di g(x)
#in un intorno della soluzione: il metodo genera una successione di iterati convergenti alla radice alfa
# ed appartenenti a questo intorno se |g'(x)|< 1 in un intorno della soluzione
plt.plot(xx,dg(xx))
plt.title('funzione dg(x)')
plt.show()


x1,it,xk=funzioni_zeri.iterazione(g,x0,tolx,nmax)
print('iterazioni= {:d}, soluzione={:e} \n\n'.format(it,x1))

#Calcolo l'ordine del metodo
ordine_iter= funzioni_zeri.stima_ordine(xk,it)
#Essendo il metodo con ordine di convergenza lineare, la costante asintotica di convergenza è data
#da |g'(alfa)| dove alfa è la radice.

C=abs(dg(x1))
print("Iterazione it={:d}, ordine di convergenza {:e}, Costante asintotica di convergenza {:e}".format(it,ordine_iter,C))

plt.plot(xx,xx,'k-',xx,g(xx))
plt.title("abs(g'(alfa))="+str(C))
Vx=[]
Vy=[]
for k in range(it):
    Vx.append(xk[k])
    Vy.append(xk[k])
    Vx.append(xk[k])
    Vy.append(xk[k+1])
    
Vy[0]=0
plt.plot(Vx,Vy,'r',xk,[0]*(it+1),'or-')
plt.show()

# Si osserva che a parità di ordine di convergenza, più piccola è la costante asintotica di convergenza,
#maggiore è la velocità del metodo.

```

## Funzioni


```python
import numpy as np
import math

'''
Il core Python non possiede la funzione sign.
La funzione copysign(a,b)  del modulo math restituisce un valore numerico che ha il valore assoluto di
a e segno di b.
Per avere il segno di un valore numerico b si può usare math.copysign(1,b)
che resistuisce 1 se b>0, -1 se b<0, 0 se b è zero
'''

def sign(x): return math.copysign(1, x)


#Bisezione
def bisez(fname,a,b,tol):
    eps=np.spacing(1)      # np.spacing(x) Restituisce la distanza tra x e il numero adiacente più vicino.
                           # np.spacing(1)  restituisce quindi l' eps di macchina.
    fa=fname(a)
    fb=fname(b)
    if sign(fa)==sign(fb):
       print('intervallo non corretto --Metodo non applicabile')
       return [],0,[]
    else:
        maxit=int(math.ceil(math.log((b-a)/tol)/math.log(2)))
        print('n. di passi necessari=',maxit,'\n');
        xk=[]
        it=0
        #while it<maxit and  abs(b-a)>=tol+eps*max(abs(a),abs(b)):
        while it<maxit and  abs(b-a)>=tol:
            c=a+(b-a)*0.5   #formula stabile per il calcolo del punto medio dell'intervallo
            xk.append(c) 
            it+=1
            if c==a or c==b:
                break
            fxk=fname(c)
            if fxk==0:
                break
            elif sign(fxk)==sign(fa):
                a=c
                fa=fxk
            elif sign(fxk)==sign(fb):
                b=c
                fb=fxk
            
            
        
        x=c
        
    return x,it,xk



def regula_falsi(fname,a,b,tol,nmax):
#Regula Falsi       
        eps=np.spacing(1)
        xk=[]
        fa=fname(a)
        fb=fname(b)
        if sign(fa)==sign(fb):
          print('intervallo non corretto --Metodo non applicabile')
          return [],0,[]
        else:
            it=0
            fxk=fname(a)
            while it<nmax and  abs(b-a)>=tol+eps*max(abs(a),abs(b)) and abs(fxk)>=tol :
                x1=a-fa*(b-a)/(fb-fa);
                xk.append(x1)
                it+=1
                fxk=fname(x1);
                if fxk==0:
                    break
                elif sign(fxk)==sign(fa):
                    a=x1;
                    fa=fxk;
                elif sign(fxk)==sign(fb):
                    b=x1;
                    fb=fxk;
                
                
            if it==nmax :
                print('Regula Falsi: Raggiunto numero max di iterazioni')
                
           
            
        return x1,it,xk
    

def corde(fname,fpname,x0,tolx,tolf,nmax):
 #Corde
        xk=[]
        m=fpname(x0)               #m= Coefficiente angolare della tangente in x0
        fx0=fname(x0)
        d=fx0/m
        x1=x0-d
        fx1=fname(x1)
        xk.append(x1)
        it=1
        while it<nmax and  abs(fx1)>=tolf and abs(d)>=tolx*abs(x1) :
           x0=x1
           fx0=fname(x0)
           d=fx0/m
           '''
           #x1= ascissa del punto di intersezione tra  la retta che passa per il punto
           (xi,f(xi)) e ha pendenza uguale a m  e l'asse x
           '''
           x1=x0-d  
           fx1=fname(x1)
           it=it+1
           xk.append(x1)
          
        if it==nmax:
            print('raggiunto massimo numero di iterazioni \n')
            
        
        return x1,it,xk
    


#Secanti
def secanti(fname,xm1,x0,tolx,tolf,nmax):
        xk=[]
        fxm1=fname(xm1);
        fx0=fname(x0); 
        d=fx0*(x0-xm1)/(fx0-fxm1)
        x1=x0-d;
        xk.append(x1)
        fx1=fname(x1);
        it=1
       
        while it<nmax and abs(fx1)>=tolf and abs(d)>=tolx*abs(x1):
            xm1=x0
            x0=x1
            fxm1=fname(xm1)
            fx0=fname(x0) 
            d=fx0*(x0-xm1)/(fx0-fxm1)
            x1=x0-d
            fx1=fname(x1)
            xk.append(x1);
            it=it+1;
           
       
        if it==nmax:
           print('Secanti: raggiunto massimo numero di iterazioni \n')
        
        return x1,it,xk
        
    
def newton(fname,fpname,x0,tolx,tolf,nmax):
#Newton
        xk=[]
        fx0=fname(x0)
        dfx0=fpname(x0)
        if abs(dfx0)>np.spacing(1):
            d=fx0/dfx0
            x1=x0-d
            fx1=fname(x1)
            xk.append(x1)
            it=0
           
        else:
            print('Newton:  Derivata nulla in x0 - EXIT \n')
            return [],0,[]
        
        it=1
        while it<nmax and abs(fx1)>=tolf and  abs(d)>=tolx*abs(x1):
            x0=x1
            fx0=fname(x0)
            dfx0=fpname(x0)
            if abs(dfx0)>np.spacing(1):
                d=fx0/dfx0
                x1=x0-d
                fx1=fname(x1)
                xk.append(x1)
                it=it+1
            else:
                 print('Newton: Derivata nulla in x0 - EXIT \n')
                 return x1,it,xk           
           
        if it==nmax:
            print('Newton: raggiunto massimo numero di iterazioni \n');
        
        return x1,it,xk
    
def stima_ordine(xk,iterazioni):
      p=[]

      for k in range(iterazioni-3):
         p.append(np.log(abs(xk[k+2]-xk[k+3])/abs(xk[k+1]-xk[k+2]))/np.log(abs(xk[k+1]-xk[k+2])/abs(xk[k]-xk[k+1])));
     
      ordine=p[-1]
      return ordine

#Newton Modificato
def newton_m(fname,fpname,x0,m,tolx,tolf,nmax):
        eps=np.spacing(1)     
        xk=[]
        #xk.append(x0)
        fx0=fname(x0)
        dfx0=fpname(x0)
        if abs(dfx0)>eps:
            d=fx0/dfx0
            x1=x0-m*d
            fx1=fname(x1)
            xk.append(x1)
            it=0
           
        else:
            print('Newton:  Derivata nulla in x0  \n')
            return [],0,[]
        
        it=1
        while it<nmax and abs(fx1)>=tolf and  abs(d)>=tolx*abs(x1):
            x0=x1
            fx0=fname(x0)
            dfx0=fpname(x0)
            if abs(dfx0)>eps:
                d=fx0/dfx0
                x1=x0-m*d
                fx1=fname(x1)
                xk.append(x1)
                it=it+1
            else:
                 print('Newton Mod: Derivata nulla   \n')
                 return x1,it,xk           
           
        if it==nmax:
            print('Newton Mod: raggiunto massimo numero di iterazioni \n');
        
        return x1,it,xk


     
def iterazione(gname,x0,tolx,nmax):
        
        xk=[]
        xk.append(x0)
        x1=gname(x0)
        d=x1-x0
        xk.append(x1)
        it=1
        while it<nmax and  abs(d)>=tolx*abs(x1) :
            x0=x1
            x1=gname(x0)
            d=x1-x0
            it=it+1
            xk.append(x1)
           
    
        if it==nmax:
            print('Raggiunto numero max di iterazioni \n')
        
        return x1, it,xk
```

# Sistemi Lineari e QR

## esame


```python
# -*- coding: utf-8 -*-
"""
Created on Tue Dec  7 12:09:22 2021

@author: HWRUser
"""

import numpy as np
import Sistemi_lineari as sl

# CHOLESKY --> autovalori > 0 e simmetrica
# LU --> determinante delle sottomatrici != 0

A = np.matrix('10 -4 4 0; -4 10 0 2; 4 0 10 2; 0 2 2 0')
print(np.linalg.eigvals(A))
#A ha autovalore negativo, quindi non ammette cholesky

B = np.matrix('5 -2 2 0; -2 5 0 1; 2 0 5 1; 0 1 1 5')
print(np.linalg.eigvals(B))
#B ha autovalori positivi e simmetrica, ammette cholesky

for i in range(1, 5):
    if np.linalg.det(A[0:i, 0:i]) == 0:
        print('Errore: determinante di A = 0')
        
#A ha tutte le sottomatrici con determinante diverso da 0, ammette LU

for i in range(1, 5):
    if np.linalg.det(B[0:i, 0:i]) == 0:
        print('Errore: determinante di B = 0')
        
#B ha tutte le sottomatrici con determinante diverso da 0, ammette LU

Ainv = np.linalg.inv(A)
Binv = np.linalg.inv(B)

P, L, U, flag = sl.LU_nopivot(A)
print(P)
print(L)
print(U)

detA = np.prod(np.diag(U))
detInvA = 1/detA

P, L, U, flag = sl.LU_nopivot(B)
print(P)
print(L)
print(U)

detB = np.prod(np.diag(U))
detInvB = 1/detB

# PA = LU --> det(P) * det(A) = det(L) * det(U)
```

## Funzioni


```python
import numpy as np

def Lsolve(L,b):
    """  
    Risoluzione con procedura forward di Lx=b con L triangolare inferiore  
     Input: L matrice triangolare inferiore
            b termine noto
    Output: x: soluzione del sistema lineare
            flag=  0, se sono soddisfatti i test di applicabilità
                   1, se non sono soddisfatti
    """
#test dimensione
    m,n=L.shape
    flag=0;
    if n != m:
        print('errore: matrice non quadrata')
        flag=1
        x=[]
        return x, flag
    
     # Test singolarita'
    if np.all(np.diag(L)) != True:
         print('el. diag. nullo - matrice triangolare inferiore')
         x=[]
         flag=1
         return x, flag
    # Preallocazione vettore soluzione
    x=np.zeros((n,1))
    
    for i in range(n):
         s=np.dot(L[i,:i],x[:i]) #scalare=vettore riga * vettore colonna
         x[i]=(b[i]-s)/L[i,i]
      
     
    return x,flag


def Usolve(U,b):
    
    """
    Risoluzione con procedura backward di Rx=b con R triangolare superiore  
     Input: U matrice triangolare superiore
            b termine noto
    Output: x: soluzione del sistema lineare
            flag=  0, se sono soddisfatti i test di applicabilità
                   1, se non sono soddisfatti
    
    """ 
#test dimensione
    m,n=U.shape
    flag=0;
    if n != m:
        print('errore: matrice non quadrata')
        flag=1
        x=[]
        return x, flag
    
     # Test singolarita'
    if np.all(np.diag(U)) != True:
         print('el. diag. nullo - matrice triangolare superiore')
         x=[]
         flag=1
         return x, flag
    # Preallocazione vettore soluzione
    x=np.zeros((n,1))
    
    for i in range(n-1,-1,-1):
         s=np.dot(U[i,i+1:n],x[i+1:n]) #scalare=vettore riga * vettore colonna
         x[i]=(b[i]-s)/U[i,i]
      
     
    return x,flag


def LUsolve(L,U,P,b):
     """
     Risoluzione a partire da PA =LU assegnata
     """
     Pb=np.dot(P,b)
     y,flag=Lsolve(L,Pb)
     if flag == 0:
         x,flag=Usolve(U,y)
     else:
        return [],flag

     return x,flag
 
    
def LU_nopivot(A):
    """
    % Fattorizzazione PA=LU senza pivot   versione vettorizzata
    In output:
    L matrice triangolare inferiore
    U matrice triangolare superiore
    P matrice identità
    tali che  LU=PA=A
    """
    # Test dimensione
    m,n=A.shape
   
    flag=0;
    if n!=m:
      print("Matrice non quadrata")
      L,U,P,flag=[],[],[],1 
      return P,L,U,flag
  
    P=np.eye(n);
    U=A.copy();
 # Fattorizzazione
    for k in range(n-1):
       #Test pivot 
          if U[k,k]==0:
            print('elemento diagonale nullo')
            L,U,P,flag=[],[],[],1 
            return P,L,U,flag

  #     Eliminazione gaussiana
          U[k+1:n,k]=U[k+1:n,k]/U[k,k]                                   # Memorizza i moltiplicatori	  
          U[k+1:n,k+1:n]=U[k+1:n,k+1:n]-np.outer(U[k+1:n,k],U[k,k+1:n])  # Eliminazione gaussiana sulla matrice
     
  
    L=np.tril(U,-1)+np.eye(n)  # Estrae i moltiplicatori 
    U=np.triu(U)           # Estrae la parte triangolare superiore + diagonale
    return P,L,U,flag

def LU_nopivotv(A):
    """
    % Fattorizzazione PA=LU senza pivot   versione vettorizzata intermedia
    In output:
    L matrice triangolare inferiore
    U matrice triangolare superiore
    P matrice identità
    tali che  LU=PA=A
    """
    # Test dimensione
    m,n=A.shape
   
    flag=0;
    if n!=m:
      print("Matrice non quadrata")
      L,U,P,flag=[],[],[],1 
      return P,L,U,flag
  
    P=np.eye(n);
    U=A.copy();
 # Fattorizzazione
    for k in range(n-1):
       #Test pivot 
          if U[k,k]==0:
            print('elemento diagonale nullo')
            L,U,P,flag=[],[],[],1 
            return P,L,U,flag

  #     Eliminazione gaussiana
          for i in range(k+1,n):
             U[i,k]=U[i,k]/U[k,k]                                   # Memorizza i moltiplicatori	  
             U[i,k+1:n]=U[i,k+1:n]-U[i,k]*U[k,k+1:n]  # Eliminazione gaussiana sulla matrice
     
  
    L=np.tril(U,-1)+np.eye(n)  # Estrae i moltiplicatori 
    U=np.triu(U)           # Estrae la parte triangolare superiore + diagonale
    return P,L,U,flag

def LU_nopivotb(A):
    """
    % Fattorizzazione PA=LU senza pivot  versione base
    In output:
    L matrice triangolare inferiore
    U matrice triangolare superiore
    P matrice identità
    tali che  LU=PA=A
    """
    # Test dimensione
    m,n=A.shape
    flag=0;
    if n!=m:
      print("Matrice non quadrata")
      L,U,P,flag=[],[],[],1 
      return P,L,U,flag
  
    P=np.eye(n);
    U=A.copy();
 # Fattorizzazione
    for k in range(n-1):
         #Test pivot 
         
         
          if U[k,k]==0:
            print('elemento diagonale nullo')
            L,U,P,flag=[],[],[],1 
            return P,L,U,flag

  #     Eliminazione gaussiana
          for i in range(k+1,n):
                U[i,k]=U[i,k]/U[k,k]
                for j in range(k+1,n):                                 # Memorizza i moltiplicatori	  
                  U[i,j]=U[i,j]-U[i,k]*U[k,j]  # Eliminazione gaussiana sulla matrice
     
  
    L=np.tril(U,-1)+np.eye(n)  # Estrae i moltiplicatori 
    U=np.triu(U)           # Estrae la parte triangolare superiore + diagonale
    return P,L,U,flag

def swapRows(A,k,p):
    A[[k,p],:] = A[[p,k],:]
    
    
def LU_pivot(A):
    """
    % Fattorizzazione PA=LU con pivot 
    In output:
    L matrice triangolare inferiore
    U matrice triangolare superiore
    P matrice di permutazione
    tali che  PA=LU
    """
    # Test dimensione
    m,n=A.shape
    flag=0;
    if n!=m:
      print("Matrice non quadrata")
      L,U,P,flag=[],[],[],1 
      return P,L,U,flag
  
    P=np.eye(n);
    U=A.copy();
 # Fattorizzazione
    for k in range(n-1):
       #Scambio di righe nella matrice U e corrispondente scambio nella matrice di permutazione per
       # tenere traccia degli scambi avvenuti
       
       #Fissata la colonna k-esima calcolo l'indice di riga p a cui appartiene l'elemento di modulo massimo a partire dalla riga k-esima
          p = np.argmax(abs(U[k:n,k])) + k
          if p != k:
              swapRows(P,k,p)
              swapRows(U,k,p)

  #     Eliminazione gaussiana
          U[k+1:n,k]=U[k+1:n,k]/U[k,k]                                   # Memorizza i moltiplicatori	  
          U[k+1:n,k+1:n]=U[k+1:n,k+1:n]-np.outer(U[k+1:n,k],U[k,k+1:n])  # Eliminazione gaussiana sulla matrice
     
  
    L=np.tril(U,-1)+np.eye(n)  # Estrae i moltiplicatori 
    U=np.triu(U)           # Estrae la parte triangolare superiore + diagonale
    return P,L,U,flag



def solve_nsis(A,B):
  # Test dimensione  
    m,n=A.shape
    flag=0;
    if n!=m:
      print("Matrice non quadrata")
       
      return
    
    Y= np.zeros((n,n))
    X= np.zeros((n,n))
    P,L,U,flag= LU_nopivot(A)
    
    if flag==0:
        for i in range(n):
            y,flag=Lsolve(L,np.dot(P,B[:,i]))
            Y[:,i]=y.squeeze(1)
            x,flag= Usolve(U,Y[:,i])
            X[:,i]=x.squeeze(1)
    else:
        print("Elemento diagonale nullo")
        X=[]
    return X    
    
```

## Es4


```python
import numpy as np
import funzioni_Sistemi_lineari as fSl
import scipy.linalg as spl
import matplotlib.pyplot as plt

xesatta=np.array([[2],[2]])

err_rel_nopivot=[]
err_rel_pivot=[]
for k in range(2,19,2):
    A= np.array([[10.0**(-k),1],[1,1]])
    b=np.array([[2+10.0**(-k)],[4]])

    P,L,U,flag, = fSl.LU_nopivot(A)

    if flag==0:
        x_nopivot,flag=fSl.LUsolve(L,U,P,b)
    else:
        print("Sistema non risolubile senza strategia pivotale")
    
    err_rel_nopivot.append(np.linalg.norm(x_nopivot-xesatta,1)/np.linalg.norm(xesatta,1))
    
    P_pivot,Lpivot,Upivot,flagpivot, = fSl.LU_pivot(A)

    if flagpivot==0:
        x_pivot,flag=fSl.LUsolve(Lpivot,Upivot,P_pivot,b)
    else:
        print("Sistema non risolubile con strategia pivotale")
        
    err_rel_pivot.append(np.linalg.norm(x_pivot-xesatta,1)/np.linalg.norm(xesatta,1))
    

plt.semilogy(range(2,19,2),err_rel_nopivot,range(2,19,2),err_rel_pivot)
plt.legend(['No pivot','Pivot'])
plt.show()
```

## esame qr


```python
"""
Created on Sat May  1 11:23:27 2021

@author: damia
"""
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as spl

#Funzioni necessarie
def Usolve(U,b):
    
    """
    Risoluzione con procedura backward di Rx=b con R triangolare superiore  
     Input: U matrice triangolare superiore
            b termine noto
    Output: x: soluzione del sistema lineare
            flag=  0, se sono soddisfatti i test di applicabilità
                   1, se non sono soddisfatti
    
    """ 
#test dimensione
    m,n=U.shape
    flag=0;
    if n != m:
        print('errore: matrice non quadrata')
        flag=1
        x=[]
        return x, flag
    
     # Test singolarita'
    if np.all(np.diag(U)) != True:
         print('el. diag. nullo - matrice triangolare superiore')
         x=[]
         flag=1
         return x, flag
    # Preallocazione vettore soluzione
    x=np.zeros((n,1))
    
    for i in range(n-1,-1,-1):
         s=np.dot(U[i,i+1:n],x[i+1:n]) #scalare=vettore riga * vettore colonna
         x[i]=(b[i]-s)/U[i,i]
      
     
    return x,flag


def metodoQR(x,y,n):
    """
    INPUT
    x vettore colonna con le ascisse dei punti
    y vettore colonna con le ordinate dei punti 
    n grado del polinomio approssimante
    OUTPUT
     a vettore colonna contenente i coefficienti incogniti
 """
 
    H=np.vander(x,n+1)
    Q,R=spl.qr(H)
    y1=np.dot(Q.T,y)
    a,flag=Usolve(R[0:n+1,:],y1[0:n+1])
    return  a



#------------------------------------------------------
#Script
#----------------------------------------------------
m=12

x=np.linspace(1900,2010,12)
y=np.array([76.0,92.0,106.0,123.0,132.0,151.0,179.0,203.0,226.0,249.0,281.0,305.0])


xmin=np.min(x)
xmax=np.max(x)
xval=np.linspace(xmin,xmax,100)

for n in range(1,4):
    a=metodoQR(x,y,n)
    residuo=np.linalg.norm(y-np.polyval(a,x))**2
    print("Norma del residuo al quadrato",residuo)
    p=np.polyval(a,xval)
    plt.plot(xval,p)
    

plt.legend(['n=1','n=2','n=3'])
 
plt.plot(x,y,'o')
```

## Metodo QR


```python
import numpy as np
import scipy.linalg as spl
from funzioni_Sistemi_lineari import Usolve
 
def metodoQR(x,y,n):
    """
    INPUT
    x vettore colonna con le ascisse dei punti
    y vettore colonna con le ordinate dei punti 
    n grado del polinomio approssimante
    OUTPUT
     a vettore colonna contenente i coefficienti incogniti
 """
 
    H=np.vander(x,n+1)
    Q,R=spl.qr(H)
    y1=np.dot(Q.T,y)
    a,flag=Usolve(R[0:n+1,:],y1[0:n+1])
    return  a

```


```python
import numpy as np
import matplotlib.pyplot as plt
from funzioni_Approssimazione_MQ import metodoQR



x = np.array([0.0004, 0.2507, 0.5008, 2.0007, 8.0013])
y = np.array([0.0007, 0.0162,0.0288, 0.0309, 0.0310]);


#Calcolo della retta di regessione
a=metodoQR(x,y,1)
residuo=np.linalg.norm(y-np.polyval(a,x))**2
print("Norma al quadrato del residuo Retta di regressione",residuo)
xmin=np.min(x)
xmax=np.max(x)
xval=np.linspace(xmin,xmax,100)
p=np.polyval(a,xval)
plt.plot(xval,p,'r-',x,y,'o')
plt.legend(['Retta di regressione', 'Dati'])
plt.show()

#Calcolo della parabola di approssimazione nel senso dei minimi quuadrati
a=metodoQR(x,y,2)
residuo=np.linalg.norm(y-np.polyval(a,x))**2
print("Norma al quadrato del residuo Polinomio di approssimazione di grado 2",residuo)
xmin=np.min(x)
xmax=np.max(x)
xval=np.linspace(xmin,xmax,100)
p=np.polyval(a,xval)
plt.plot(xval,p,'r-',x,y,'o')
plt.legend(['Polinomio di approssimazione di grado 2', 'Dati'])
plt.show()

```

# Interpolazione

## Funzione


```python
import numpy as np
def plagr(xnodi,k):
    """
    Restituisce i coefficienti del k-esimo pol di
    Lagrange associato ai punti del vettore xnodi
    """
    xzeri=np.zeros_like(xnodi)
    n=xnodi.size
    if k==0:
       xzeri=xnodi[1:n]
    else:
       xzeri=np.append(xnodi[0:k],xnodi[k+1:n])
    
    num=np.poly(xzeri) 
    den=np.polyval(num,xnodi[k])
    
    p=num/den
    
    return p



def InterpL(x, f, xx):
     """"
        %funzione che determina in un insieme di punti il valore del polinomio
        %interpolante ottenuto dalla formula di Lagrange.
        % DATI INPUT
        %  x  vettore con i nodi dell'interpolazione
        %  f  vettore con i valori dei nodi 
        %  xx vettore con i punti in cui si vuole calcolare il polinomio
        % DATI OUTPUT
        %  y vettore contenente i valori assunti dal polinomio interpolante
        %
     """
     n=x.size
     m=xx.size
     L=np.zeros((n,m))
     for k in range(n):
        p=plagr(x,k)
        L[k,:]=np.polyval(p,xx)
    
    
     return np.dot(f,L)
```

## Es


```python

import numpy as np
from funzioni_Interpolazione_Polinomiale import InterpL
import matplotlib.pyplot as plt
#nodi del problema di interpolazione 

T=np.array([-55, -45, -35, -25, -15, -5, 5, 15, 25, 35, 45, 55, 65])
L=np.array([3.7, 3.7,3.52,3.27, 3.2, 3.15, 3.15, 3.25, 3.47, 3.52, 3.65, 3.67, 3.52])

# punti di valutazione per l'interpolante
xx=np.linspace(np.min(T),np.max(T),200);
pol=InterpL(T,L,xx);

pol42=InterpL(T,L,np.array([42]))
pol_42=InterpL(T,L,np.array([-42]))

plt.plot(xx,pol,'b--',T,L,'r*',42,pol42,'og',-42,pol_42,'og');
plt.legend(['interpolante di Lagrange','punti di interpolazione','stima 1', 'stima2']);
plt.show()

```

# Integrazione

## esame


```python
"""
Created on Sat Jun  5 08:29:40 2021

@author: damia
"""
import numpy as np
import matplotlib.pyplot as plt
import math

#funzioni per l'integrazione Simposon Composita e con ricerca automatica del numero di N
#di sottointervalli

def SimpComp(fname,a,b,n):
    h=(b-a)/(2*n)
    nodi=np.arange(a,b+h,h)
    f=fname(nodi)
    I=(f[0]+2*np.sum(f[2:2*n:2])+4*np.sum(f[1:2*n:2])+f[2*n])*h/3
    return I


def simptoll(fun,a,b,tol):

    Nmax=4096
    err=1
    
    N=1;
    IN=SimpComp(fun,a,b,N);
    
    while N<=Nmax and err>tol :
        N=2*N
        I2N=SimpComp(fun,a,b,N)
        err=abs(IN-I2N)/15
        IN=I2N
 
    
    if N>Nmax:
        print('Raggiunto nmax di intervalli con simptoll')
        N=0
        IN=[]
 
    return IN,N

#Funzioni per l'interpolazione di Lagrange
def plagr(xnodi,k):
    """
    Restituisce i coefficienti del k-esimo pol di
    Lagrange associato ai punti del vettore xnodi
    """
    xzeri=np.zeros_like(xnodi)
    n=xnodi.size
    if k==0:
       xzeri=xnodi[1:n]
    else:
       xzeri=np.append(xnodi[0:k],xnodi[k+1:n])
    
    num=np.poly(xzeri) 
    den=np.polyval(num,xnodi[k])
    
    p=num/den
    
    return p



def InterpL(x, f, xx):
     """"
        %funzione che determina in un insieme di punti il valore del polinomio
        %interpolante ottenuto dalla formula di Lagrange.
        % DATI INPUT
        %  x  vettore con i nodi dell'interpolazione
        %  f  vettore con i valori dei nodi 
        %  xx vettore con i punti in cui si vuole calcolare il polinomio
        % DATI OUTPUT
        %  y vettore contenente i valori assunti dal polinomio interpolante
        %
     """
     n=x.size
     m=xx.size
     L=np.zeros((n,m))
     for k in range(n):
        p=plagr(x,k)
        L[k,:]=np.polyval(p,xx)
    
    
     return np.dot(f,L)
 
# Script principale
tol=1e-08
x=np.zeros((6,))
y=np.zeros((6,))
N=np.zeros((6,))
fig=1

#Funzione integranda
f= lambda x: 2/math.pi*(5.5*(1-np.exp(-0.05*x))*np.sin(x**2))


for i in range(0,6):
    #Estremo destro dell'intervallo di integrazione
    x[i]=0.5+2*i
   
    #Costruisco punti nell'intervallo [0,x[i]] in cui valutare e poi disegnare la funzione integranda
    xi=np.linspace(0,x[i],100)
    plt.subplot(2,3,fig)
    plt.plot(xi,f(xi))
    plt.legend([ 'x['+str(i)+']='+str(x[i])])
    fig+=1
    #Calcolo il valore dell'integrale i-esimo  con a=0 e b= x[i]con la precisione richiesta
    y[i],N[i]=simptoll(f,0,x[i],tol)
 
plt.show()
xx=np.linspace(min(x),max(x),100)

#Calcolo il polinomio che interpola le coppie (x,y)
pol=InterpL(x, y, xx)
plt.plot(xx,pol,x,y,'ro')
plt.legend(['Polinomio interpolante','Nodi di interpolazione'])
plt.show()
print("Numero di sottointervalli per ciascuni il calcolo di ciascun integrale \n",N)

```


```python


import numpy as np
import f_Interpol_Polinomial as fip
import f_Integrazi_numerica as fin
import matplotlib.pyplot as plt

A = 1
B = 3
n = 3

f = lambda x: x-np.sqrt(x-1)

punti = np.linspace(A, B, 50)
nodi = np.linspace(A, B, n+1)
pol = fip.InterpL(nodi, f(nodi), punti)
print(pol)

def fp(val):
    return fip.InterpL(nodi, f(nodi), val)

plt.plot(punti, f(punti), 'r', punti, pol, 'b', nodi, f(nodi), 'go')
plt.show()

tol = 1e-5
I1 = 2.114381916835873
I2 = 2.168048769926493

I1a, N1 = fin.simptoll(f, A, B, tol)
I2a, N2 = fin.simptoll(fp, A, B, tol)

print(N1)
print(N2)
print(np.abs(I1a - I1))
print(np.abs(I2a - I2))
```


```python

import numpy as np
import f_Interpol_Polinomial as ip
import matplotlib.pyplot as plt

f = lambda x: 1/(1+900*x**2)

a = -1
b = 1

punti = np.linspace(a, b, 50)

for n in range(5, 35, 5):
    x = np.array(list())
    for i in range(1, n+2):
        x = np.append(x, -1+(2*(i-1)/n))
            
    pe = ip.InterpL(x, f(x), punti)
    re = np.abs(f(punti) - pe)
    plt.plot(punti, re)
    plt.show()
    
for n in range(5, 35, 5):
    x = np.array(list())
    for i in range(n+1, 0, -1):
        x = np.append(x, np.cos(((2*i-1)*np.pi)/(2*(n+1))))
            
    pe = ip.InterpL(x, f(x), punti)
    re = np.abs(f(punti) - pe)
    plt.plot(punti, re)
    plt.show()
```


```python

import numpy as np
import f_Interpol_Polinomial as fip
import matplotlib.pyplot as plt

f = lambda x: np.cos(np.pi * x) + np.sin(np.pi * x)
x = np.array([1.0, 1.5, 1.75])

a = 0
b = 2
xx = np.linspace(a, b, 50)

p = fip.InterpL(x, f(x), xx)

plt.plot(xx, p, 'r', xx, f(xx), 'g', x, f(x), 'o')
plt.legend(['Polinomio', 'Funzione', 'Nodi'])
plt.show()

xo = np.array([0.75])
r = np.abs(f(0.75) - fip.InterpL(x, f(x), xo))
print('Resto: ', r)

x = np.array([0.75, 1.0, 1.5, 1.75])
p2 = fip.InterpL(x, f(x), xx)

plt.plot(xx, p2, 'r', xx, f(xx), 'g', x, f(x), 'o')
plt.legend(['Polinomio 2', 'Funzione', 'Nodi'])
plt.show()
```

## Funzioni


```python
import numpy as np

def TrapComp(fname,a,b,n):
    h=(b-a)/n
    nodi=np.arange(a,b+h,h)
    f=fname(nodi)
    I=(f[0]+2*np.sum(f[1:n])+f[n])*h/2
    return I

    
def SimpComp(fname,a,b,n):
    h=(b-a)/(2*n)
    nodi=np.arange(a,b+h,h)
    f=fname(nodi)
    I=(f[0]+2*np.sum(f[2:2*n:2])+4*np.sum(f[1:2*n:2])+f[2*n])*h/3
    return I


def traptoll(fun,a,b,tol):

    Nmax=2048
    err=1
    
    N=1;
    IN=TrapComp(fun,a,b,N);
    
    while N<=Nmax and err>tol :
        N=2*N
        I2N=TrapComp(fun,a,b,N)
        err=abs(IN-I2N)/3
        IN=I2N
 
    
    if N>Nmax:
        print('Raggiunto nmax di intervalli con traptoll')
        N=0
        IN=[]
 
    return IN,N


def simptoll(fun,a,b,tol):

    Nmax=2048
    err=1
    
    N=1;
    IN=SimpComp(fun,a,b,N);
    
    while N<=Nmax and err>tol :
        N=2*N
        I2N=SimpComp(fun,a,b,N)
        err=abs(IN-I2N)/15
        IN=I2N
 
    
    if N>Nmax:
        print('Raggiunto nmax di intervalli con traptoll')
        N=0
        IN=[]
 
    return IN,N
```

## Es1


```python
import sympy as sym
import Funzioni_Integrazione as FI 
from sympy.utilities.lambdify import lambdify
import numpy as np
import matplotlib.pyplot as plt
scelta=input("Scegli funzione ")

x=sym.symbols('x')
scelta_funzione = {
        '1': [x**10,0.0,1.0],
        '2': [sym.asin(x),0.0,1.0],
        '3': [sym.log(1+x), 0.0,1.0]
}
fx,a,b=scelta_funzione.get(scelta)

Iesatto=float(sym.integrate(fx,(x,a,b)))


f= lambdify(x,fx,np)
N=[1, 2, 4, 8, 16, 32 ,64 ,128, 256]
i=0
InT=[]
InS=[]
for n in N:
    InT.append(FI.TrapComp(f,a,b,n))
    InS.append(FI.SimpComp(f,a,b,n))
  
ET=np.zeros((9,))
 
ES=np.zeros((9,))


ET=np.abs(np.array(InT)-Iesatto)/abs(Iesatto)

ES=np.abs(np.array(InS)-Iesatto)/abs(Iesatto)

plt.semilogy(N,ET,'ro-',N,ES,'b*-')
plt.legend(['Errore Trapezi Composita', 'Errore Simpson Composita'])
plt.show()
```

## Es2 toll


```python
import sympy as sym
import Funzioni_Integrazione as FI 
from sympy.utilities.lambdify import lambdify
import numpy as np
import matplotlib.pyplot as plt
scelta=input("Scegli funzione ")

x=sym.symbols('x')
scelta_funzione = {
        '1': [sym.log(x),1.0,2.0],
        '2': [sym.sqrt(x),0.0,1.0],
        '3': [sym.Abs(x), -1.0,1.0]
}
fx,a,b=scelta_funzione.get(scelta)

Iesatto=float(sym.integrate(fx,(x,a,b)))


f= lambdify(x,fx,np)


tol=1e-6

IT,NT=FI.traptoll(f,a,b,tol)

print("Il valore dell'integrale esatto e' ", Iesatto)
if NT>0:
    print("Valore con Trapezi Composito Automatica ",IT," numero di suddivisoini ",NT)

IS,NS=FI.simptoll(f,a,b,tol)

if NS>0:    
    print("Valore con Simpson Composito Automatica ",IS," numero di suddivisoini ",NS)
```

c871f7a49375fc5e64f5733cf1b51721ef3e8f9b
ipynb
Jupyter Notebook
Untitled11.ipynb
sergiorolnic/github-slideshow
589404d7a9f217135d19ef661ce203851a02b6b5
MIT
Untitled11.ipynb
sergiorolnic/github-slideshow
589404d7a9f217135d19ef661ce203851a02b6b5
MIT
2020-10-10T21:39:27.000Z
2020-10-11T09:55:57.000Z
Untitled11.ipynb
sergiorolnic/github-slideshow
589404d7a9f217135d19ef661ce203851a02b6b5
MIT
Qwen/Qwen-72B
1. YES
2. YES