I want to solve, in MatLab, a linear system (corresponding to a PDE system of two equations written in finite difference scheme). The action of the system matrix (corresponding to one of the diffusive terms of the PDE system) reads, symbolically (
u is one of the unknown fields,
n is the time step,
j is the grid point):
The above matrix has to be intended as
A, where A*U^n+1 = B is the system.
U contains the 'u' and the 'v' (second unknown field of the PDE system) alternatively: U = [u_1,v_1,u_2,v_2,...,u_J,v_J]. So far I have been filling this matrix using
diag in the following expensive way:
E=zeros(2*J,1); E(1:2:2*J) = 1; E(2:2:2*J) = 0; Dvec=zeros(2*J,1); for i=3:2:2*J-3 Dvec(i)=D_11((i+1)/2); end for i=4:2:2*J-2 Dvec(i)=D_21(i/2); end A = diag(Dvec)*spdiags([-E,-E,2*E,2*E,-E,-E],[-3,-2,-1,0,1,2],2*J,2*J)/(dx^2);`
and for the solution
[L,U]=lu(A); y = L\B; U(:) =U\y;
B is the right hand side vector.
This is obviously unreasonably expensive because it needs to build a JxJ matrix, do a JxJ matrix multiplication, etc.
Then comes my question: is there a way to solve the system without passing MatLab a matrix, e.g., by passing the vector
Dvec or alternatively directly
D_22? This would spare me a lot of memory and processing time!