Risoluzione tramite fem del problema di propagazione di calore in un corpo in Python 3.0

Ricostruiamo il puzzle completo dato che tutti i pezzi che lo compongono li abbiamo già sul tavolo, avendoli definiti negli articoli precedenti.

Il problema

Consideriamo il seguente problema differenziale espresso correttamente e posto in forma forte. Ciò significa che il problema è posto sotto forma di un sistema composto oltre che dalla equazione differenziale stessa, anche dalle condizioni iniziali e al contorno – di Dirichelet – necessarie per poter arrivarvi a determinare una soluzione univoca. Il problema è lo stesso dell’esercizio 6.8 della 2a edizione del testo Partial Differential Equations di Mark S. Gockenbach.

$$\left\{\begin{matrix}
\rho c \frac{\partial u}{\partial t} -\kappa \frac{\partial^2 u}{\partial x^2}=f(x,t) , 0<x<100, t>0
\\ u(x,0)=\psi(x), 0<x<100 \\u(0,t)=0, t>0
\\ u(100,t)=0, t>0 \end{matrix}\right.$$

Esso descrive la distribuzione del calore interna ad una barra di ferro di densità di materia $\rho=7.88 [\frac{kg}{m^3}]$, calore specifico (o capacità termica per unità di massa) $c=0.437 [\frac{J}{K kg}]$ e conducibilità termica $\kappa = 0.836[\frac{W}{m K}]$. La barra ha lunghezza $l=100$ $cm$ ed è scaldata internamente ad un tasso $f(x,t)=10^{-8}tx(100-x)^2 [\frac{W}{cm^3}]$ mentre entrambi gli estremi sono posti ad una temperatura costante (bagno termico) di $0 °C$. Per semplicità è possibile fattorizzare le tre costanti all’interno della diffusività termica $\alpha=\frac{\kappa}{\rho c}$

L’esercizio proposto rappresenta un tipico problema di diffusione. La forma differenziale impiegata è di tipo parabolico poiché la sua equazione caratteristica ha determinante nullo. Prima di proseguire, è simpatico riportare quelli che sono i casi notevoli:

• l’equazione di Poisson nel caso di sorgenti costanti $\bigtriangledown u=\rho$
• l’equazione di Laplace per sorgenti nulle $\bigtriangledown u=0$

L’obiettivo che l’esercizio ci chiede di perseguire è quello di trovare la distribuzione di temperatura su tutta la barra trascorsi 180 secondi di tempo dall’istante iniziale. Per farlo, utilizzo tutti gli step elencati al dettaglio negli articoli precedenti i quali, in seguito, verranno argomentati e implementati nelle sezioni di codice che li compete. Il codice nella sua interezza verrà allegato a fine pagina in formato .py. Proseguiamo.

Teorema di Lax-Milgram

Cominciamo anticipando che, un problema ellittico su un dominio $\Omega \in R^d$ ha esistenza e unicità della soluzione assicurate dal teorema di Max-Milgram se e solo se :

• il dominio V è un opportuno spazio di Hilbert
• a è una forma bilineare, continua e coerciva
• f è un funzionale lineare e continuo da V in $R$

Librerie, costanti e condizioni iniziali

Il primo passo da compiere è quello di includere le librerie necessarie per plottare e maneggiare numericamente le funzioni. Scelgo di plottare usando la nonna di tutte le librerie grafiche matplotlib.pyplot che, dopo 10 anni, rimane ancora molto valida e completa. Includo la libreria numpy per poterne usare funzioni e metodi e per manipolare vettori e matrici, e la libreria scipy per integrare e derivare con facilità.

#librerie 
import matplotlib.pyplot as plt     
import numpy as np     
import scipy as sp   

Qui apro e chiudo una piccola parentesi: Scipy, acronimo di scientific python è una libreria davvero utile, e di recente è stata sottoposta ad un notevole miglioramento. Dalla versione 1.4 vi è stata implementata la funzione quad_vec che permette di integrare funzioni vettoriali. Attualmente l’ultima versione è la 1.6.2 uscita il 24/03/21. Per conoscerne la versione è sufficiente eseguire il seguente codice:

"""
import scipy as sp
from scipy import version
version=sp.version.version
print("version scipy=",version)
"""

Successivamente dichiaro tutte le costanti che caratterizzano il sistema $\rho$, $\kappa$ e $c$, e gli estremi del dominio spazio-temporale di definizione $t_0=0$, $t_{fin}=180$, $x_0=0$ e $x_{fin}=l=100$

#costanti
ro=7.88
c=0.437
k=0.836
A=k/ro*c      #coefficiente di diffusività termica []

#estremi del dominio spaziale e temporale
x0=0
l=100
t0=0
tfin=180

Definisco anche le funzione initialTemperatureDistro ovvero la distribuzione di temperatura iniziale interna alla barra (zero gradi celsius ovunque) e la funzione rhs ovvero la forzante f(x,t) al problema considerato, il nostro termine di sorgente.

#distribuzione iniziale della temperatura
def initialTemperatureDistro(x): 
    return 0

#right hand side della equazione: è il calore immesso sul sistema 
def rhs(x,t):
    return 10**(-8)*t*x*(100-x)**2

Forma debole

Come abbiamo visto negli articoli precedenti (ad esempio vedi qui), risolvere numericamente un problema differenziale impiegando il metodo degli elementi finiti richiede una previa trasformazione del problema in forma debole. Per farlo, moltiplico l’equazione differenziale per una funzione di test $v(x)$ ed integro. Semplifico andando ad elidere alcuni termini per via delle condizioni al contorno imposte al sistema e, a questo punto, non voglio più risolvere il problema iniziale ma il seguente: $$\int_{0}^{100} \left \{ \rho c \frac{\partial u(x,t)}{\partial t}v(x) + \kappa \frac{\partial u(x,t)}{\partial x}\frac{\partial v(x)}{\partial x} \right \}dx=\int_{0}^{100}f(x,t)v(x)dx, \forall v \in V$$

Il metodo Galerkin

Un problema posto in forma debole ammette l’applicazione del metodo Galerkin (vedi qui) per il quale si approssima lo spazio infinito-dimensionale $S$ di partenza con uno spazio finito-dimensionale $S_n$ dato da $$S_n=\left \{ p:[0,l]\rightarrow \mathbb{R}, p(0)=p(l)=0 : \mathrm{p\ sia\ continua\ e\ lineare\ a\ tratti} \right \}$$ a sua volta generato da funzioni di base opportune $$S_n=span\left \{ \phi_1,\phi_2,…,\phi_{n-1} \right \}$$

Triangolazione del dominio

Partiamo partizionando il dominio di definizione, cioè l’intervallo dell’asse reale che va da 0 a 100 cm, in sotto-intervalli di ampiezza costante. Chiamo $n=100$ il numero di sottodomini, $h=100/n$ l’ampiezza del dominio, $x_i=ih$ la posizione dell’elemento finito e $i=0,1,…,n$ l’indice che identifica l’elemento finito

#partizione del dominio spaziale
nstep=100
h=1.0*(l/nstep)
xarray=np.linspace(x0,l,nstep+1)
dx=xarray[1]-xarray[0]

#stampa su terminale
print("dx= ",dx)
print("xarray= ",xarray)

Discretizzo, poi, anche il dominio temporale. In questo caso, il singolo quanto temporale è $dt=0.6$.

#discretizzazione del dominio temporale
ntstep=300
dt=(tfin-t0)/ntstep
tarray=np.linspace(t0,tfin,ntstep+1)
dt=tarray[1]-tarray[0]

#stampa su terminale
print("dt= ",dt)
print("tarray= ",tarray)

Base di polinomi a pezzi

Ora è giunto il momento di definire la base utilizzata per approssimare la funzione in ogni sottospazio $S_n$ di $S$. Essendo il numero di nodi per sottospazio pari a $n=2$, la base scelta all’interno di ognuno di essi, sarà lineare (a tratti) e sarà la base standard fatta di di triangoli (roof functions) $\left \{ \phi_1,\phi_2,…,\phi_{n-1} \right\}$. Matematicamente si deve avere

$$\phi_i(x)=\left\{\begin{matrix}
\frac{1}{h}[x-(i-1)h] & x_{i-1}<x<x_i\\
-\frac{1}{h}[x-(i+1)h]& x_i<x<x_{i+1}\\
0& altrimenti
\end{matrix}\right.$$

la cui derivata diventa

$$\frac{\mathrm{d} \phi_i}{\mathrm{d} x}(x)=\left\{\begin{matrix}
\frac{1}{h} & x_{i-1}<x<x_i\\
-\frac{1}{h}& x_i<x<x_{i+1}\\
0& altrimenti
\end{matrix}\right.$$

All’interno del codice definisco preventivamente una funzione ausiliaria diff(j,x) che valuta il rapporto tra la differenza fra la x considerata e l’elemento j-esimo della base (cui compete la posizione j*h sull’asse delle x) e l’ampiezza d’intervallo h. Questo passaggio risulterà certamente superfluo, io l’ho aggiunto per evitare un possibile refuso matematico nella funzione myPiecewise cioè la funzione che va a definire le roof functions vere e proprie

#funzione ausiliaria per la base di polinomi
def diff (j,x):
    return (x-(j)*(h))/h 

#base di polinomi lineari a pezzi 
def myPiecewise(x,i):
    p=np.piecewise(x, [x<(i-1)*h,((i-1)*h<=x)&(x<(i*h)),(i*h<=x)&(x<(i+1)*h),x>=(i+1)* h], [0,diff(i-1,x),1-diff(i,x),0])
    return p

Questo codice ricrea la base seguente: una base fatta da 100 triangoli, uno per ogni sotto-intervallo in cui si è divisa la lunghezza della barra di 100cm. Ciascun triangolo funge da base nel rispettivo sottoinsieme in cui è definito. Ciò significa che il suo valore sarà zero ovunque tranne che all’interno di quello stesso intervallo e, in particolare, varrà zero sui due nodi agli estremi e varrà $h=1$ al vertice, nel punto in cui la funzione è discontinua.

Per comprendere meglio quel che sta succedendo, di seguito vengono riportati 10 elementi nel caso in cui lo l’intervallo di definizione sia suddiviso diviso in soli 10 sotto-intervalli. Come si può vedere ad occhio nudo, i nodi sull’asse delle x sono 10+1 mentre i triangoli sono solo 9. Questo per evitare che un triangolo, cioè una base rimanga incompleto agli estremi. Inoltre, fra ogni coppia di nodi, la funzione è lineare. il suo gradiente è costante. La tipica forma a triangolo o a tetto appare solo quando si accostano tra loro non 2 ma 3 nodi.

Di seguito si considera un solo elemento, per evidenziare il fatto che ciascuna base a forma triangolare assume valore diverso da zero solo all’interno del proprio sottospazio.

Invece del solito plot statico, per divertimento ho deciso di usare plt.ion() per realizzare un grafico dinamico. La base di polinomi viene disegnata ad uno ad uno all’interno del ciclo for.

plt.ion()
fig=plt.figure()
#plt.figure()

a=[]
for k in range(nstep-1):
    a.insert(k,vmyPiece(k+1,xarray))
    plt.plot(xarray,a[k]) 
    plt.xlabel('x')
    plt.ylabel('a')
    plt.grid(True)
    plt.pause(0.1)
plt.waitforbuttonpress(0)
plt.close()

In generale è sempre comodo vettorizzare le funzioni per fare in modo che qualsiasi trasformazione che voglio svolgere su dei vettori sia fatta nella loro interezza, sul vettore globale e non su ogni singolo elemento che lo compone. Così facendo si risparmiano risorse utili. Per farlo utilizzo la funzione vectorize definita all’interno della libreria numpy

#vettorizzazione
vmyPiece=np.vectorize(myPiecewise)

Proiezioni ortogonali

Per ogni tempo t la funzione $u(.,t)$ sarà opportunatamente approssimata all’interno di ogni sottospazio $S_n$ e potremo riscriverne la soluzione cercata come somma del prodotto delle funzioni di base $\phi_i(x)$ per degli opportuni coefficienti $\alpha_i$ che dipendono esclusivamente da $t$: $$u_n(x,t)=\sum_{i=1}^{n-1}\alpha_i(t)\phi_i(x)$$

Ri-esprimo il problema andando, anzitutto, a sostituire la funzione incognita che appare come integranda con la sua approssimazione data da una somma finita di termini $u_n(x,t)$ $$\int_{0}^{100} \left \{ \rho c \frac{\partial u_n(x,t)}{\partial t}v(x) + \kappa \frac{\partial u_n(x,t)}{\partial x}\frac{\partial v(x)}{\partial x} \right \}dx=\int_{0}^{100}f(x,t)v(x)dx, \forall v \in S_n$$

dopodiché, essendo le $v(x)$ funzioni qualsiasi del sottospazio $V$, le vado a sostituire con una base $\phi_j(x)$ per $V$

$$\int_{0}^{100} \left \{ \rho c \frac{\partial u_n(x,t)}{\partial t}\phi_i(x) + \kappa \frac{\partial u_n(x,t)}{\partial x}\frac{\partial \phi_i(x)}{\partial x} \right \}dx=$$ $$\int_{0}^{100}f(x,t)\phi_i(x)dx, i=1,2,…,n-1$$

Sostituendo $u_n(x,t)$ con la sua forma esplicita abbiamo

$$\sum_{j=1}^{n-1}\frac{\mathrm{d} \alpha_j}{\mathrm{d} t}(t)\int_{0}^{100} \rho c \phi_j(x)\phi_i(x)dx + \sum_{j=1}^{n-1}\alpha_j(t)\int_{0}^{100}\kappa \frac{\partial \phi_j(x)}{\partial x}\frac{\partial \phi_i(x)}{\partial x}dx=$$ $$\int_{0}^{100}f(x,t)\phi_i(x)dx, i=1,2,…,n-1$$

E’ questa l’equazione algebrica cui vogliamo arrivare! Equazione che sarà risolta coi metodi numerici classici. Ma prima, esprimiamola in una forma che ci è più familiare.

Il termine di destra è il prodotto della forzante per la base di polinomi, integrati su tutto il dominio.

$$f(t)=\begin{bmatrix}
\int_{0}^{l}f(x,t)\phi_1(x)dx\\
\int_{0}^{l}f(x,t)\phi_2(x)dx\\
.\\
.\\
\int_{0}^{l}f(x,t)\phi_{n-1}(x)dx\
\end{bmatrix}$$

dove $$f_i(t)=\int_{0}^{100}f(x,t)\phi_i(x)dx=t\frac{h^2(60000i-200h-1200i^2h+6h^2i^3+3ih^2)}{6\times 10^8},$$ $$i=1,2,…,n-1$$

Lo esprimo definendo la funzione scalare myRhs come

def myRhs(x,t,i):
    return myPiecewise(x,i)*rhs(x,t)

e introducendo la funzione fvalues come quel il vettore che ha per elementi gli integrali di myRhs su tutto il dominio. Per integrare utilizzo la funzione quad interna alla libreria scipy. fvalues opera così: ogni volta che viene chiamata prende come argomento la variabile scalare t che esprime l’istante temporale considerato, chiama la funzione quad da scipy.integrate e definisce al suo interno un array vuoto sigma di lunghezza nstep-1 (un vettore). Per ogni indice i che assume tutti i valori interi a partire da 0 e fino a 100, all’elemento i di sigma viene assegnato il valore corrispondente all’integrale definito da 0 a 100 di myRhs in dx. fvalues infine restituisce l’array di primitive.

def fvalues(t):
    from scipy.integrate import quad
    sigma=np.empty([nstep-1])
    for i in range (nstep-1):
        sigma[i]=quad(myRhs,0,l,args=(t,i+1))[0]
    return sigma

N.B. un array di tipo np.array parte dall’elemento numero zero quindi dire che ha lunghezza nstep-1 (e dimensione (1,nstep-1)) dove nstep-1=100-1=99 significa assegnarli come elemento anche lo zero e quindi avrà comunque lunghezza 100

N.B.B. in quad args=(,) sta ad indicare parametri extra che non dipendono dalla variabile di integrazione. Porre [0] al termine restituisce errore nullo.

Matrici tridiagonali

Nei termini di sinistra è possibile identificare la matrice di massa

$$M_{ij}=\int_{0}^{l}\rho c \phi_j(x)\phi_i(x)dx$$

e la matrice di rigidezza

$$K_{ij}=\int_{0}^{l}\kappa \frac{\mathrm{d} \phi_j}{\mathrm{d} x}(x)\frac{\mathrm{d} \phi_i}{\mathrm{d} x}(x)dx$$

Per implementarle nel codice utilizzo una loro fondamentale proprietà: il fatto di essere sparse, cioè il fatto che la maggior parte dei loro elementi risulti essere nulla ed, in particolare, il fatto che gli elementi non nulli risultino essere posti lungo la diagonale principale e le due diagonali secondarie. Questo fatto è molto importante. L’aver scelto una base triangolare ha reso queste due matrici tridiagonali e simmetriche. Ciò significa che invertirle non richiederà un grande sforzo per il calcolatore (anche se l’inversa non sarà più tridiagonale). Nel codice, prima definisco la matrice generale tridiagMatrix che ricostruisce una matrice tridiagonale nel seguente modo: prende come argomenti e restituisce

#definizione della matrice tridiagonale
def tridiagMatrix(T,x,y,z,k1=-1, k2=0, k3=1):
    a = [x]*(T-abs(k1)); b = [y]*(T-abs(k2)); c = [z]*(T-abs(k3))
    return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3)

Di seguito creo le due matrici considerando che:

• per i termini sulla diagonale principale $$K_{ii}\int_{0}^{100}\kappa(\frac{\mathrm{d} \phi_i}{\mathrm{d} x}(x))^2dx=\frac{2\kappa}{h}, i=1,2,…,n-1$$ $$M_{ii}\int_{0}^{100}\rho c (\phi_i(x))^2dx=\frac{2h \rho c }{3}, i=1,2,…,n-1$$
• per i termini sulle diagonali secondarie $$K_{i,i+1}=\int_{0}^{100}\kappa \frac{\mathrm{d} \phi_i}{\mathrm{d} x}(x)\frac{\mathrm{d} \phi_{i+1}}{\mathrm{d} x}(x)dx=-\frac{\kappa}{h}, i=1,2,…,n-1$$ $$M_{i,i+1}=\int_{0}^{100}\rho c \phi_i(x) \phi_{i+1}(x)dx=\frac{h \rho c }{6}, i=1,2,…,n-2$$

#creazione della matrice di rigidezza (stiffness matrix)
K=tridiagMatrix(nstep-1,-0.836,2*0.836,-0.836)                  
print("K=",K,type(K),np.shape(K)) 

#creazione della matrice di massa (mass matrix)
M=tridiagMatrix(nstep-1,(h*ro*c)/6,(2*h*ro*c)/3,(h*ro*c)/6)      
print("M=",M,type(M),np.shape(M)) 

Ed ecco le matrici stampate su terminale

Equazioni ai coefficienti

Infine, chiamando $\alpha(t)$ i coefficienti incogniti

$$\alpha(t)=\begin{bmatrix}
\alpha_1(t)\\
\alpha_2(t)\\
.\\
.\\
\alpha_{n-1}(t)
\end{bmatrix}$$

non resta che risolvere l’equazione algebrica che segue:

$$M\frac{\mathrm{d} \alpha}{\mathrm{d} t}+K\alpha=f(t)$$

Da un problema differenziale al contorno ai valori iniziali espresso in forma forte siamo arrivati al dover risolvere un problema di Cauchy: una equazione differenziale ordinaria ai valori iniziali

$$M\frac{\mathrm{d} \alpha}{\mathrm{d} t}+K\alpha=f(t), t>t_0$$

$$\alpha(t_0)=\alpha_0$$

con $$\alpha_0=\begin{bmatrix}
\psi(x_1)\\
\psi(x_2)\\
.\\
.\\
\psi(x_{n-1})
\end{bmatrix}$$

In generale, la tecnica appena usata, di risolvere una PDE dipendente dal tempo andando ad integrare l’equivalente sistema semidiscreto di ODE viene chiamata tecnica delle linee.

Risoluzione della equazione algebrica associata

Come risolvere questo problema? Una volta arrivati alla equazione ai coefficienti il più è fatto, si possono usare diverse tecniche risolutive a questo punto. Proviamo ad esprimere il problema in forma discreta:

$$\alpha^{i+1}=\alpha^i + \Delta t M^{-1}(-K\alpha^i+f(t_i))$$

Poiché è sconveniente invertire una matrice tridiagonale ($M^{-1}$ diventerebbe densa) conviene optare per una via alternativa. Poniamoci il problema di risolvere l’equaizone equivalente $$\alpha^{i+1}=\alpha^i+\Delta t s^i$$ con $$Ms^i=-K\alpha^i+f(t_i)$$

Un metodo conveniente da utilizzare è quello di invocare la funzione solve_banded di scipy, una funzione ritagliata appositamente per risolvere le espressioni algebriche che coinvolgono matrici tridiagonali. Riportiamo la documentazione inerente la funzione presa direttamente dal sito : https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html

scipy.linalg.solve_banded
#Solve the equation a x = b for x, assuming a is banded matrix.

Prima di essere chiamata solve_banded necessita di 3 argomenti obbligatori, di cui uno deve essere prima definito, e sono:

• l and u – rappresentano il numero di diagonali inferiori e superiori della matrice a bande in cui appaiono elementi non nulli
• ab – una nuova matrice che costruisco a partire dalle diagonali non nulle della matrice iniziale. Poiché nella espressione algebrica da risolvere compaiono due matrici, ab sarà una combinazione di due nuove matrici costruite a partire da M e K. Nel seguito vedremo come, nel frattempo ci basta sapere che ab deve avere dimensione pari a [nr delle diagonali non nulle]x[numero degli elementi della diagonale principale]
• b – la righthand side della espressione

Per costruire la matrice ab ho prima ricreato separatamente 6 array vuoti di lunghezza pari ad una colonna della matrice M (la stessa di K) che poi ho intenzione di popolare con gli elementi delle diagonali diverse da zero.

am=np.empty(len(M))
ak=np.empty(len(K))
um=np.empty(len(M))
uk=np.empty(len(K))
lm=np.empty(len(M))
lk=np.empty(len(K))

Per popolare gli array utilizzo un ciclo for effettuato sull’indice che rappresenta il singolo elemento interno ad ogni colonna, o riga, di M (K) che è una matrice quadrata. Gli elementi di M (K) presi rispettivamente dalla diagonale superiore, dalla diagonale principale e da quella inferiore vengono riassegnati agli array vuoti

• am (ak) rappresenta l’array che ha per componenti gli elementi sulla diagonale principale di M (K)
• um (uk) è l’array che ha per componenti gli elementi sulla diagonale secondaria superiore di M (K)
• lm (lk) è l’array che ha per componenti gli elementi sulla diagonale secondaria inferiore di M (K).

for i in range (len(M)):
    for j in range (len(M[i])):
        if (i==j):
            am[j]=M[i][j]
            ak[j]=K[i][j]
            #print("M[i][j]=",M[i][j])
        if (i==j-1):
            um[j]=M[i][j]
            uk[j]=K[i][j]
            #print("M[i][j]=",M[i][j])
        if (i==j+1):
            lm[j]=M[i][j]
            lk[j]=K[i][j]
            #print("M[i][j]=",M[i][j])

Infine cucio assieme gli array creati ad hoc grazie alla funzione vstack di numpy e stampo su terminale il risultato ottenuto, giusto perché non fa mai male avere una verifica diretta di ogni step.

mb=np.vstack((um,am,lm))
kb=np.vstack((uk,ak,lk))
print("mb=",mb,type(mb),np.shape(mb))
print("kb=",kb,type(kb),np.shape(kb))

vstack incolla in una unica matrice gli array che le vengono passati come argomento, rispettandone l’ordine.

Backward Euler

Invoco il metodo di Eulero all’indietro, un metodo implicito (a differenza del metodo di Eulero classico) perché la funzione incognita cambia andamento molto rapidamente – è una funzione stiff.

Il fatto che sia un metodo implicito significa che in ogni step di valutazione del valore assunto della funzione in un punto (un nodo della griglia), il uso valore dipende da tutte le altre variabili in quel punto e non solamente dal valore delle stesse in punti precedenti. Quindi queste dovranno essere rivalutate step per step.

mback=mb+dt*kb
alpha=np.zeros(nstep-1)

Finalmente, invoco solve_banded e lo faccio assegnandone il valore alla variabile temporanea coeff.

for t in range(0,ntstep+1):
    N=-np.dot(K,alpha)
    coeff=solve_banded((1,1),mback,-np.dot(K,alpha)+fvalues(t*dt))
    alpha=alpha+dt*coeff

Definisco la soluzione ricercata come una funzione ad un sol parametro in ingresso:

def sol(x):
    sol=0
    for i in range(nstep-1):
        sol=sol+myPiecewise(x,i+1)*alpha[i]
    return sol

Vettorizzo la funzione soluzione

vsol=np.vectorize(sol)

E la chiamo assegnando alla variabile s il valore che la funzione assume in tutti i punti della griglia spaziale grid. Averla prima vettorizzata permette di fare in modo che per ogni punto della griglia l’intera la funzione venga valutata complessivamente.

grid=np.linspace(0,100,101)
s=vsol(grid)

Plot

plt.plot(grid,s)
plt.xlabel('x')
plt.ylabel('a')
plt.grid(True)
plt.axis((0,100,0,4)) 
plt.show()   
#plt.pause(0.1)    
plt.waitforbuttonpress(0)
plt.close()