Appunti Programmazione - calcolo numerico - Python

# HELP DELLA FUNZIONE THOMAS

# INPUT:

# D vettore con gli elementi della diagonale principale della matrice dei coefficienti

# S vettore con gli elementi sulla diagonale superiore della matrice dei coefficienti

# A vettore con gli elementi sulla diagonale inferiore della matrice dei coefficienti

# B vettore dei termini noti

# OUTPUT:

# X vettore della soluzione

#

def Thomas(D,S,A,B):

# import librerie

from numpy import zeros, size

# inizializzazione array

N = size(B)

V= zeros(N-1)

U = zeros(N)

alpha = zeros(N-1)

Y = zeros(N)

X = zeros(N)

# Fattorizzazione LU

V = S

U[0] = D[0]

for i in range(N-1):

alpha[i] = A[i]/U[i]

U[i+1] = D[i+1]-alpha[i]*V[i]

# soluzione del sistema lineare

Y[0] = B[0]

for i in range(N-1):

Y[i+1] = B[i+1]-alpha[i]*Y[i]

X[N-1] = Y[N-1]/U[N-1]

for i in range(N-2,-1,-1):

X[i]=(Y[i]-V[i]*X[i+1])/U[i]

# output

return X

# input

import numpy as np

D=np.array([1,2,3,4])

S=np.array([1,2,3])

A=np.array([1,2,3])

B=np.array([1,2,3,4])

# soluzione:

soluzione=Thomas(D,S,A,B)

print soluzione

# METODO UPWIND

# Soluzione dell'equazione del trasporto con il METODO UPWIND

#

# Input:

# c: velocità di trasporto

# T0, T: intervallo temporale

# a, b: intervallo di discretizzazione spaziale

# M: numero di passi nel tempo

# N: numero di nodi (spaziali) interni

#

# Output:

# Xi, Tj: nodi

# U[0:N+1,0:M]: soluzione approssimata

#

# Function:

# phi_t0(x): condizione iniziale

# psi_xl(t): condizione al bordo sinistro

# psi_xr(t): condizione al bordo destro

import numpy as np

import sys

import matplotlib.pyplot as plt

# Condizione iniziale

phi_t0 = lambda x: np.sin(x)

# Condizione al bordo sinistro

psi_xl = lambda t: -np.sin(2.0+4.0*t)

# Condizione al bordo destro

psi_xr = lambda t: np.sin(2.0-4.0*t)

# DATI DI INPUT:

# Lettura dati di input (come stringa)

dati_input = raw_input("Inserisci i valori c T0 T a b M N separati da uno spazio: ")

# Trasforma la stringa di input in numeri

dati = []

for valore in dati_input.split():

try: dati.append(float(valore))

except: print("I numeri devono essere separati da uno spazio!!")

sys.exit()

if len(dati) != 7:

print("Devi inserire 7 numeri!!")

sys.exit()

# Dati di input del problema

c, T0, T, a, b, M, N = dati[0], dati[1], dati[2], dati[3], dati[4], int(dati[5]), int(dati[6])

print "c: ", c, "T0: ", T0, "T: ", T, "a: ", a, "b: ", b

print "M: ", M, "N: ", N

# Passo temporale

k = (T-T0)/M

print('Passo temporale: ', k)

# Passo spaziale

h = (b-a)/(N+1)

print('Passo spaziale: ', h)

# Nodi

Tj = np.linspace(T0,T,M+1)

Xi = np.linspace(a,b,N+2)

# Numero di Courant

alpha = c*k/h

print('Numero di Courant: ', alpha)

# Inizializzazione del vettore delle approssimazioni

U = np.zeros((N+2,M+1))

# Condizione iniziale

for i in xrange(N+2):

U[i,0] = phi_t0(Xi[i])

# Condizione al bordo sinistro

U[0,:] = psi_xl(Tj)

# Condizione al bordo destro

U[N+1,:] = psi_xr(Tj)

# Algoritmo del metodo upwind

for j in xrange(M):

U[1:N+1,j+1] = (1-alpha)*U[1:N+1,j] + alpha*U[0:N,j]

# oppure faccio con due cicli for:

# Algoritmo del metodo upwind

#for j in xrange(M):

# for i in xrange(1,N+1):

# U[i,j+1] = (1-alpha)*U[i,j] + alpha*U[i-1,j]

# stessa identica cosa

# Grafico 2D della soluzione

plt.figure(1)

plt.clf()

plt.plot(Xi,U)

plt.xlabel('x')

plt.ylabel('u(x,t)')

plt.show()

# Libreria per i grafici 3D

from mpl_toolkits.mplot3d import Axes3D

# Griglia

Tgrid, Xgrid = np.meshgrid(Tj, Xi)

fig = plt.figure(2)

plt.clf()

ax = Axes3D(fig)

ax.plot_surface(Xgrid,Tgrid,U)

ax.set_xlabel('x')

ax.set_ylabel('t')

plt.show()

# Esercizio risolto con metodo upwind

def UpWind(c, T0, k, M, a, b, N):

from numpy import zeros

# Passo spaziale

h = (b-a)/(N+1)

# Numero di Courant

alpha = c*k/h

# Nodi

Xi=zeros(N+2)

Tj=zeros(M+1)

for i in xrange(N+2):

Xi[i]=a+i*h

for j in xrange(M+1):

Tj[j]=j*k

# Inizializzazione del vettore delle approssimazioni

U = zeros((N+2,M+1))

# Condizione iniziale

for i in xrange(N+2):

U[i,0] =phi_t0(Xi[i])

# Condizione al bordo sinistro

for j in xrange(M+1):

U[-1,j] = psi_xl(Tj[j])

# Algoritmo del metodo upwind

for j in xrange(M):

U[1:N+1,j+1] = (1-alpha)*U[1:N+1,j] + alpha*U[0:N,j]

print "Xi:", Xi, '\n' "Tj:",Tj, '\n' "Uij", U

return Xi, Tj, U

import matplotlib.pyplot as plt

import numpy as np

from math import pi

# DATI DI INPUT:

k=0.1

a=-1.

b=1.

c=1./3.

T0=0.

M=6

N=8

# FUNZIONI ESTERNE:

def phi_t0(x):

return abs(np.sin(pi*x))

def psi_xl(t):

return abs(np.sin(pi*t/3))

# richiamo funzione:

a=UpWind(c, T0, k, M, a, b, N)

Xi=a[0]

Tj=a[1]

U=a[2]

# GRAFICI:

# Grafico 2D della soluzione

plt.figure(1)

plt.clf()

plt.plot(Xi,U)

plt.xlabel('x')

plt.ylabel('u(x,t)')

plt.show()

# Libreria per i grafici 3D

from mpl_toolkits.mplot3d import Axes3D

# Griglia

Tgrid, Xgrid = np.meshgrid(Tj, Xi)

fig = plt.figure(2)

plt.clf()

ax = Axes3D(fig)

ax.plot_surface(Xgrid,Tgrid,U)

ax.set_xlabel('x')

ax.set_ylabel('t')

plt.show()

# cicli di iterazione

# metodo iterativo di jacobi

# Soluzione di un sistema lineare con il metodo di Jacobi

#

#

# Input:

# A[1:N,1:N], B[1:N]: matrice dei coefficienti e termine noto

# X0[1:N]: approssimazione iniziale

# epsilon: accuratezza richiesta

# Niter_max: numero massimo di iterazioni

# Output:

# Xnew[1:N]: nuova approssimazione

# Errk[1:Niter_max]: norma della differenza tra due approssimazioni successive

# Niter: numero di iterazioni necessarie per raggiungere l’accuratezza richiesta

#

def Jacobi(A,B,X0,epsilon,Niter_max):

import numpy as np

from linalg import norm

# inizializzazione

N = np.size(B)

Errk = np.zeros(Niter_max)

CJ= np.zeros((N,N))

QJ= np.zeros((N,N))

# metodo jacobi

CJ=-np.dot(np.diag(1./np.diag(A)),(np.triu(A)+np.tril(A)))

QJ=np.dot(np.diag(1./np.diag(A)),B)

# ciclo di iterazione

# prima devo dare delle condizioni da cui iniziare su err e icount. se le condizioni richieste

dal ciclo while sono verificate entra nel ciclo while e procede

err = 10*epsilon # voglio entrare nel ciclo while quindi parto da una condizione per cui

riesco a entrare. come per esempio 10*epsilon, che è sicuramente maggiore di epsilon (avrei potuto

mettere una condizione qualsiasi)

icount = 1 # idem di err: parto da una condizione arbitraria tale che io possa entrare nel ciclo

while

# X0 è definito tramite i dati di input della funzione. entra nel ciclo con quell'X0.

while ( (err > epsilon) and (icount < Niter_max) ): # se le condizioni su err e icount sono

verificate entra nel ciclo while

Xnew = np.dot(CJ,X0)+QJ # definisce l'approssimazione nuova mediante la

precedente, che si chiama X0. alla prima iterazione sarà l'X0 dei dati di input

err = norm(Xnew-X0) # calcola l'errore a questa iterazione (appross nuova - appross

vecchia) X0= np.copy(Xnew) # butta via l'approssimazione vecchia e la sostituisce con quella

nuova (in questo caso la vecchia si chiama X0. alla prima iterazione ho l'X0 dei dati input, dalla

seconda in poi diventa via via l'approssimazione al passo precedente)

Errk[icount-1] = err # memorizza l'errore all'iesima iterazione

icount += 1 # butta via il vecchio indice di iterazione icount e ne crea uno nuovo ari a

icount+1 (avrei potuto scrivere anche icount=icount +1)

# a questo punto riprova le condizioni del while sul nuovo icount e se sono verificate

rinizia da capo il ciclo while con i nuovi dati

# ovviamente si blocca anche se arriva la condizione su err

return Xnew, Errk, icount-1

# Cj=-D^-1(L+U) --> non usare l'inversa !!!!!!

# Qj=D^-1 B

# D=diag(diag(A)) --> il primo diag dice che deve creare una matrice diagonale, il secondo diag

dice deve estrarre gli elementi dellla diagonale di A (ricordarsi l'import delle librerie)

# U=triu(A)

# L=tril(A)

# D^-1= diag(1./diag(A)) --> l'inversa di una matrice o vettore in python si può fare facendo

1/matrice.

# attenzione all'ordine perchè se prima creo la matrice diagonale con elementi diagonali pari alla

diagonale di A e poi faccio l'inversa 1/matrice mi ritrovo a dividere per zero in molti elementi

# Soluzione di un sistema lineare con il metodo iterativo di Jacobi

#

#

# Input:

# A[1:N,1:N], B[1:N]: matrice dei coefficienti e termine noto

# X0[1:N]: approssimazione iniziale

# epsilon: accuratezza richiesta

# Niter_max: numero massimo di iterazioni

# Output:

# Xnew[1:N]: nuova approssimazione

# Errk[1:Niter_max]: norma della differenza tra due approssimazioni successive

# Niter: numero di iterazioni necessarie per raggiungere l’accuratezza richiesta

#

def Jacobi(A,B,X0,epsilon,Niter_max):

import numpy as np

from numpy.linalg import norm

# inizializzazione

N = np.size(B)

Errk = np.zeros(Niter_max)

CJ= np.zeros((N,N))

QJ= np.zeros((N,N))

# metodo jacobi

CJ=-np.dot(np.diag(1./np.diag(A)),(np.triu(A)+np.tril(A)))

QJ=np.dot(np.diag(1./np.diag(A)),B)

# ciclo di iterazione

err = 10*epsilon

icount = 1

Xold=X0

while ( (err > epsilon) and (icount < Niter_max) ):

Xnew = np.dot(CJ,Xold)+QJ

err = norm(Xnew-Xold)

Xold= np.copy(Xnew)

Errk[icount-1] = err

icount += 1

return Xnew, Errk, icount-1

import numpy as np

import sys

# INPUT

dati_input = raw_input("Inserisci i valori di N, epsilon e Niter_max separati da uno spazio: ")

dati = []

for valore in dati_input.split():

try: dati.append(float(valore))

except: print("I numeri devono essere separati da uno spazio!!")

sys.exit()

N=int(dati[0])

epsilon=dati[1]

Niter_max=int(dati[2])

print "N:", N, "epsilon:", epsilon, "Niter_max:", Niter_max

# vettore dei termini noti:

inputB = raw_input("Inserisci gli elementi del vettore dei termini noti separati da uno spazio: ")

elementiB=[]

for valore in inputB.split():

try: elementiB.append(float(valore))

except: print("I numeri devono essere separati da uno spazio!!")

sys.exit()

B=np.array(elementiB)

# Matrice dei coefficienti A:

inputelementiA = raw_input("Inserisci gli elementi della matrice dei coe