I try to solve a heat diffusion problem on tetrahedral finite elements with nodal heat sources, which depend on the solution vector, in MATLAB. The nonlinear equation system looks like this:
BU' + AU = q(T)
with B being the heat capactiy matrix, A being the conductivity matrix, q being the source terms and U being the Temperature. I use an Adams-Bashforth/Trapezoid Rule predictor-corrector scheme with a Picard iteration followed by a time step control. The temperature for the source terms is evaluated exactly between the last time step's temperature and the predictor's temperature. Here is a simplified version of the predictor-corrector code. The calculation of the sources is a function.
% predictor
K0 = t(n)-t(n-1);
Upre(dirichlet) = u_d_t(coordinates(dirichlet,:));
Upre(FreeNodes) = U(FreeNodes,n) + (dt/2)*((2+dt/K0)*U_dt(FreeNodes,3) - (dt/K0)*U_dt(FreeNodes,1)); % predictor step
Uguess = Upre; % save as initial guess for Picard iteration
% corrector with picard iteration
while res >= picard_tolerance
T_theta = Uguess*theta + (1-theta)*U(:,n);
b = q(T_theta);
% Building right-hand side vector (without Dirichlet boundary conditions yet)
rhs = ((2/dt)*B*U(:,n) + B*U_dt(:,1))+b;
% Applying Dirichlet Boundary Conditions to the Solution Vector
Ucor(dirichlet) = u_d_t(coordinates(dirichlet,:));
rhs = rhs - ((2/dt)*B+A)*Ucor;
% Solving the linearized system using the backslash operator
% P*U(n+1) = f(Un) => U(n+1) = P\f(Un)
Ucor(FreeNodes) = ((2/dt)*B(FreeNodes,FreeNodes)+A(FreeNodes,FreeNodes))\rhs(FreeNodes);
res = norm(Uguess-Ucor);
Uguess = Ucor;
U(:,n+1) = Ucor;
end
As you can see I use the backslash operator to solve the system. The non-linearity of the system should not be to bad. However, with increasing time steps the picard method converges more slowly and eventually stops to converge altogether. I need much bigger time steps though, so I put the whole corrector step into a function and tried to solve it with fsolve instead to see if I achieve quicker convergence. Unfortunately fsolve seems to never even finish the first time step. I suppose I did not configure the options for fsolve correctly. Can anyone tell me, how to configure fsolve for large sparse nonlinear systems ( We are talking about thousands upt to hundredthousands of equations). Or is there maybe a better solution than fsolve for this problem? Help and - as I am not an expert or computational engineer - explicit advice is greatly appreciated!
