Mathematical Engineering - Matematica Numerica

Exercise - 13 - Solution

Divided by topic

MATEMATICA NUMERICA A.A. 2018 - 2019 Ingegneria Matematica Prof. A. Quarteroni Prof. A. Manzoni, Dr. I. Fumagalli Esercitazione 13 - Soluzione Problemi nonlineari  Interpolazione polinomiale Esercizio 1 Si consideri il seguente problema ai limiti, per un dato valore di0: (u00 +u3 =f(x); x2(0;1); u(0) = 0; u(1) = 0:(1) Siaf(x) = 2 +x3 (1x)3 . Con tale denizione dif(x), la soluzione esatta del problema risulta uex( x) =x(1x). Supponiamo di voler approssimare tale soluzione con il metodo degli elementi niti. La formulazione debole del problema (1) risulta data da : trovareu2H1 0(0 ;1)tale che a(u; v) +Z 1 0u 3 v=F(v)8v2V dovea(u; v) =Z 1 0u 0 v0 dx; F(v) =Z 1 0f vdx: Si consideri una partizione dell'intervallo(0;1)inN= 100elementi, si introduca lo spazioV h= X1 h\ V degli elementi niti lineari, indicando conf' igN i=1una sua base associata ai nodi 0 =x 0< x 1< : : : < x N< xN+1= 1 . Il problema discretizzato risulta : trovareu h2 V htale che a(u h; v h) + Z 1 0u 3 hv h= F(v h) 8v h2 V h: Algebricamente, tale problema è equivalente al seguente sistema diNequazioni non lineari inNincognite : Au+g(u) =f(2) doveu= (u 1; : : : ; u N)T indica il vettore delle incognite (valori nodali diu h( x) =P N j=1u j' j( x)), e per i; j= 1; : : : ; N, Aij= a(' j; ' i) ; f i= F(' i) ; g i( u) =Z 1 0u 3 h' idx: In maniera equivalente, il problema (2) può essere riscritto comeG(u) =Au+g(u)f=0: Siamo interessati a risolvere il sistema non lineare (2) con il seguente metodo di punto sso : Au( k+1) =fg(u( k) ): In modo equivalente,u( k+1) = (u( k) ) =A 1 (fg(u( k) )):1 a)Prendendo spunto dalle precedenti esercitazioni, costruire un data le per il problema lineare (u00 =f(x); x2(0;1); u(0) = 0; u(1) = 0:(3) Si considericome unico parametro variabile. b)Partendo dai template forniti, completare il main leEser13_DiffNonlinearReact_FixedPoint e la funzione[u, FE_SPACE, MESH, DATA] = Elliptic_nonlinear_1D_Solver_fixed_point_Eser13 (FE_SPACE, MESH, DATA, u0, param) aggiungendo le linee di codice richieste per implementare il metodo di punto sso proposto. Tale funzione prende in ingresso le strutture dati e la soluzione inizialeu0ottenute da una precedente chiamata al solutore lineareElliptic1D_Solver, oltre al vettoreparam, e restituisce le stesse strutture insieme alla soluzioneudel problema nonlineare. Per assemblare a ogni passo il termine nonlineare, serve una funzione dix,t,param, che può essere costruita approssimando(u( k) h)3 mediante il suo interpolante lineare a tratti, con la seguente sintassi : DATA.force = @(x,t,param)( interp1(MESH.vertices, u_k.^3, x) ) ; Come soluzione inizialeu0si consideri la soluzione del problema lineare (3). NB: si posizioni lo scriptElliptic_nonlinear_1D_Solver_fixed_point_Eser13.mcom- pletatonella cartella Code. c)SiaN= 100. Risolvere il problema con il metodo di punto sso proposto per= 1,= 10e= 20. Utilizzare il metodopcgper risolvere il sistema lineare a ciascuna iterazione, una tolleranza di10 6 per il metodo di punto sso e per il metodopcg. Denire il massimo numero di iterazioni per il metodo di punto sso epcgpari a100. Commentare i risultati ottenuti. Soluzione 1 Si tratta di implementare in Matlab il seguente metodo di punto sso per un sistema diNequazioni non lineari inNincognite : Au( k+1) =fg(u( k) );(4) considerando come dato inizialeu(0) la soluzione del problema lineareAu(0) =f, versione algebrica del problema (3). Come criterio di arresto per l'algoritmo, possiamo utilizzare quello basato sull'incremento, ossia fermarci quando u( k) u( k1) < tol, essendotol >0una tolleranza scelta. Per rendere più eciente l'assemblaggio del termine nonlneare, si considera un'approssimazione di(u h( x))3 mediante il suo interpolante lineare composito, come segue : (u h( x))3 =0 @N X j=1u j' j( x)1 A3 N X j=1u 3 j' j( x) =:e Uh( x): Il risultato di questa approssimazione è una funzione elementi niti, che si può utilizzare nel termine noto aggiuntivo per il sistema (4) : gi( u) =Z 1 0u 3 h' idx Z 1 0N X j=1u 3 j' j' idx =Z 1 0e Uh' idx: (5) Si noti che la denizione die Uhcome funzione di xpermette di assemblare gli elementig i( u( k) )del termine nonlineare in maniera analoga al vettoref: ricordiamo infatti che fi=Z 1 0f ' idx: Pertanto, possiamo assemblare il vettoreg(u( k) )tramiteAssembler_1D:2 DATA.force = @(x,t,param)( interp1(MESH.vertices, u_old.^3, x ) ); [~, guk, ~] = Assembler_1D(MESH, DATA, FE_SPACE); guk = guk(MESH.internal_dof, 1); a)Occorre innanzitutto risolvere il sistema lineare (ottenuto prendendo= 0nella formulazione del problema) per determinare il dato iniziale per le iterazioni di Newton. Per risolvere tale problema scriviamo usare il seguente data leEser13_DiffNonlinearReact_FixedPoint_data(in cuil'unico parametro variabile)data.flag_dirichlet = [1 2]; data.flag_neumann = []; data.flag_robin = []; data.bcDir = @(x,t,param)(0. *x);data.bcNeu = @(x,t,param)(0. *x);data.bcRob_alpha = @(x,t,param)(0. *x);data.bcRob_gamma = @(x,t,param)(0. *x);data.bcRob_fun = @(x,t,param)(0. *x);data.force = @(x,t,param)(2+ param(1) *x.^3. *(1-x).^3);data.diffusion = @(x,t,param)(1+0. *x);data.transport = @(x,t,param)(0+0. *x);data.reaction = @(x,t,param)(0+0. *x);e la seguente porzione di main le Eser13_DiffNonlinearReact_FixedPoint(non è richiesta la vi- sualizzazione della soluzioneu0, qui riportata per completezza). Si osservi che il parametroalpha passato al solutore lineare è diverso da zero, in modo da implementare il corretto termine forzante.close all clear all clc addpath(Code); %% %datafile data_file_name =Eser13_DiffNonlinearReact_FixedPoint_data; ...1)definizionedidominioenumerodielementi... a = 0; b = 1; N = 100; %Elementifinitilineari fem =P1; ...2)definizionedeiparametridelsistema... alpha = 20; ...3)costruzionedellamesh... [vertices,elements,boundaries] = mesh1D(a,b,N); ...4)calcolodellasoluzionedelproblemalineareausiliario... [u0, FE_SPACE, MESH, DATA] = Elliptic1D_Solver ...(elements,vertices,boundaries,fem,data_file_name,alpha); ...5)graficodellasoluzionedelpblinausiliario... figure(11) feplot(u0,MESH.nodes, MESH.elements,fem,r) xlabel(x); ylabel(u0(x)); grid on hold on 3 u_ex = @(x) x. *(1-x);plot(MESH.nodes,u_ex(MESH.nodes),g--); hold off legend({u^{(0)},u_{ex}}) title([\alpha=,num2str(alpha)]); b)Per calcolare la soluzione del problema non lineare, è suciente introdurre nel main le Eser13_DiffNonlinearReact_FixedPointle sole istruzioni necessarie alla chiamata di Elliptic_nonlinear_1D_Solver_fixed_point_Eser13; nel caso in esame, tracciamo anche il graco della soluzione esatta, in modo da valutare l'accuratezza della soluzione approssimata....6)soluzionedelproblemanonlineare, ...utilizzandolasoluzionedelpasso4comesoluzioneiniziale... [uh,~,MESH] = Elliptic_nonlinear_1D_Solver_fixed_point_Eser13 ...(FE_SPACE, MESH, DATA, u0, alpha); ...7)graficodellasoluzioneesatta ...ediquellanumericadelproblemanonlineare... figure(12) feplot(uh,MESH.nodes, MESH.elements,fem,r) xlabel(x); ylabel(uh(x)); grid on hold on u_ex = @(x) x. *(1-x);plot(MESH.nodes,u_ex(MESH.nodes),g--); hold off legend({u^{(k_{end})},u_{ex}}) title([\alpha=,num2str(alpha)]); Si noti come la soluzione del problema lineare, calcolata con la funzione Elliptic_1D_Solver, sia pas- sata alla funzioneElliptic_nonlinear_1D_Solver_fixed_point_Eser13come argomento di input. Quest'ultima funzione può essere implementata come segue, tenendo conto delle speciche fornite.function [u, FE_SPACE, MESH, DATA] = Elliptic_nonlinear_1D_Solver_fixed_point_Eser13...(FE_SPACE, MESH, DATA, u0, param) if nargin < 5DATA.param = []; else DATA.param = param; end %operatorilineari A = FE_SPACE.A_in; F = FE_SPACE.F_in; u_D = u0(MESH.Dirichlet_dof); tol = 1e-6; nmax = 100; uk_int = u0(MESH.internal_dof,1); uk = u0; u_old = zeros(size(MESH.vertices,2),1); nit = 1; alpha = DATA.param(1); while norm(uk-u_old)>tol...aggiornamentodellevariabili... u_old = uk; u_old_int = uk(MESH.internal_dof); ...assemblaggiodig(uk)... DATA.force = @(x,t,param)( interp1(MESH.vertices, u_old.^3, x ) ); [~, guk, ~] = Assembler_1D(MESH, DATA, FE_SPACE); 4 guk = guk(MESH.internal_dof, 1); ...iterazionedipuntofisso... %uk_int=A\(F-param(1) *guk);[uk_int,flag,~,iter] = pcg ( A, F - param(1) *guk, tol, nmax);iter if flag ~= 0fprintf(ERRORE(it.ptofisso=%d):pcgnonconvergein%diterazioni\n, nit, nmax); return end nit=nit+1; uk = zeros(MESH.numNodes,1); uk(MESH.internal_dof) = uk_int; uk(MESH.Dirichlet_dof) = u_D; if nit>nmaxdisp(ERRORE:nonconverge); break ;end end if nittolu_old_int = uk_int; u_old = uk; DATA.force = @(x,t,param)(interp1(MESH.vertices,alpha *u_old.^2,x));[~, w_uk, ~] = Assembler_1D(MESH, DATA, FE_SPACE); w_newton = w_uk(MESH.internal_dof); %Newton JF = A + diag(3 *w_newton);Fres = A *u_old_int + w_newton. *u_old_int - F;5 Figure 1  Soluzione esatta e approssimata (all'ultimo passo del metodo di punto sso) per il problema nonlineare (1).%deltau=JF(-Fres); [deltau, flag, ~, iter_pcg] = pcg(JF,-Fres,0.5,1000); if flag==0disp([ConvergenzadiPCGinnum2str(iter_pcg)iterazioni]); else disp(PCGnonharaggiuntoconvergenza->uscita) u = uk; return end norm(deltau) uk_int = u_old_int + deltau; nit=nit+1; uk = zeros(MESH.numNodes,1); uk(MESH.internal_dof) = uk_int; uk(MESH.Dirichlet_dof) = u_D; if nit>nmaxdisp( *errore *nonconverge);break ;end end if nit