Mathematical Engineering - Matematica Numerica

Exercise - 05 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 5 Metodi diretti per sistemi lineari Esercizio 1 (P)A partire dai modelli di main le e dataleEser05_AdvectionDiffusion_template.me Eser05_AdvectionDiffusion_data_template.m, implementare inPoliFEMla soluzione del seguente pro- blema di diusione e trasporto, usando elementi niti lineari : ((u0 )0 + u0 = 0;0< x 08x2Rn conx6 =0. A tale scopo si utilizza il seguente criterio : condi- zione necessaria e suciente anchè una matriceA2Rn n sia denita positiva è che le sottomatrici principaliA idi Adi ordinei= 1; : : : ; n(cioè quelle ottenute limitandoAalle primeirighe e colonne) abbiano determinante positivo (criterio diSylvester). Si verica con i seguenti comandi che la matrice Ain questione soddis le ipotesi di questo criterio ed è anche una matrice simmetrica : if isequal(A,A')disp('A e'' simmetrica') na = size(A, 1); for i = 1:naif (det(A(1:i,1:i))>0)if (i == na) disp('A e'' definita positiva') end else error('A non e'' definita positiva') end end else disp('A non e'' simmetrica') end A e' simmetrica A e' definita positiva Sarà dunque possibile ottenere la fattorizzazione LU. Consideriamo ora la matriceB; si può osservare che è una matrice a dominanza diagonale stretta per colonne, cioèjB(j; j)j>P n i=1;i6 =jj B(i; j)jper ogni colonnaj. Con i seguenti comandi verichiamo che la matriceBsoddisfa questa ipotesi per l'esistenza della fattorizzazione LU : db = diag(B); nb = size(B,1); for j = 1:nb r(j) = abs(db(j)) - sum(abs(B([1:j-1,j+1:nb], j))); end if (r > 0)disp('B e'' a dominanza diagonale stretta per colonne')5 end B e' a dominanza diagonale per colonne Consideriamo inne la matriceC. Essa non presenta alcuna caratteristica particolare, per cui veri- chiamo che la fattorizzazione LU sia applicabile mediante la condizione necessaria e suciente. Con i seguenti comandi verichiamo pertanto che i suoi minori principali di ordine da1adn1hanno determinante diverso da zero : i = 1; dc = 1; nc = size(C,1); while ( (i < nc) && (dc ~=0 ) )dc = det(C(1:i, 1:i)); i = i+1; end if (dc == 0)disp('non posso applicare la decomposizione LU ad C') else disp('posso applicare la decomposizione LU a C') end posso applicare la decomposizione LU a C b)Si vedalugauss.m. c)Calcoliamo la fattorizzazione LU delle precedenti matrici utilizzando il metodo di Gauss. MatriceA: [LA, UA] = lugauss(A) LA =1.0000 0 0 0.0200 1.0000 0 0.0600 -0.0100 1.0000 UA =50.0000 1.0000 3.00000 5.9800 -0.0600 0 0 0.8194 Verichiamo inne l'esattezza della fattorizzazione tramite i comandi A - LA*UA ans = 0 0 0 0 0 0 0 0 0 Consideriamo ora la matriceB: [LB, UB] = lugauss(A)6 LB = 1.0000 0 0 0.0600 1.0000 0 0.2000 0.1906 1.0000 UB =50.0000 1.0000 10.00000 19.9400 0.4000 0 0 67.9238 Verichiamo anche in questo caso l'esattezza della fattorizzazione tramite il comando B - LB*UB ans =0 0 0 0 0 0 0 0 0 Consideriamo Inne la matriceC: [LC, UC] = lugauss(C) LC = 1.0000 0 0 0.7143 1.0000 0 0.1429 -0.5000 1.0000 UC =7.0000 8.0000 9.00000 -1.7143 -3.4286 0 0 3.0000 e verichiamo perCl'esattezza della fattorizzazione : C - LC*UC ans =0 0 0 0 0 0 0 0 0 d)Si vedano le funzionifwsub.mebksub.m e)Costruiamo il termine notob: b = A * ones(na,1); Per risolvere il sistema lineareAx=busiamo le funzionifwsub.mebksub.mcon i seguenti comandi : sol = ones(3,1); b = A * sol;7 y = fwsub(LA, b); x = bksub(UA, y) x = 1.0000 1.0000 1.0000 f )La norma2dell'errore relativo residuo viene calcolata con il seguente comando err_rel = norm(sol - x)/norm(x) err_rel =1.4333e-16 e la norma2del residuo relativo si puo ottenere con il comando res_nor = norm(b - A*x)/norm(b) res_nor = 1.6267e-17 Esercizio 3 Un generatore di corrente genera una correnteI 0che uisce nelle nresistenze elettricheR 1; R 2; : : : R nin parallelo, come schematizzato in Fig. 2. Nota la correnteI 0, vogliamo calcolare le correnti I k, con k= 1; : : : ; n, che uiscono neglinrami. La tensioneVai capi di ogni resistenzaR ksegue la legge di Ohm: V=R kI k: Per ottenere lenequazioni che ci servono per risolvere il sistema, utilizziamo le leggi diKirchhoalle maglie e ai nodi dello schema elettrico. Per la legge di Kirchho ai nodi otteniamo l'equazione : I0= I 1+ I 2+


+I n;Figure 2  Schema elettrico del problema8