CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Main CFD Forum (https://www.cfd-online.com/Forums/main/)
-   -   Solving u_t = Au_x +Bu_y (https://www.cfd-online.com/Forums/main/92273-solving-u_t-au_x-bu_y.html)

VincentD September 8, 2011 09:51

Solving u_t = A(u_x +u_y) in MATLAB
 
I would like to solve the following equation:
u_t =A(u_x + u_y)
Where A is a matrix with size equal to u. I use MATLAB to do the computation.
I want the scheme to be explicit in time and implicit in space.

Before I adress this equation I will first make it a step more simple.
Lets say we want to solve
u_t = a(u_x +u_y)
Where a is a scalar.

Rewriting this equation gives:

u*-dt a(u*_x +u*_y) = u^n

We can write this equation as a matrix product in the following way:

(I+dt a(D2D))u* =u^n

Where I is the identity matrix and D2D is the matrix approximating -d/dx -d/dy. We would this D2D matrix in the following way.

We define the difference matrix D (approximating -d/dn, 1D)

spdiags([1 1 -1 ;ones(n-2,1)*[1 0 -1]; 1 1 -1],-1:1,n,n)'/(2*h); (Neumann BCs)

And we can create D2D by:

kron(speye(ny),D2D(nx-1))+kron(D2D(ny),speye(nx-1))

The full matrix I+dt a(D2D) is than

kron(speye(nx-1),speye(ny))+dt*a*kron(speye(ny),D(nx-1,hx,2,1))+kron(D(ny,hy,3,3),speye(nx-1))) (don't mind the arguments in D, first one denotes size)

Alright, well this works. But now we make it a little bit more difficult, namely we change a to be a matrix instead of a scalar

u_t =A(u_x + u_y)

The matrix equation becomes
(I+dt A(D2D))u* =u^n

Now we take that nx=ny=10
Now A is a matrix of 9x10 and D2D has size 90x90 and same for I (90x90).

I need to blow matrix A up in such a way that it's also square with a size (90x90), however if I agian use the kron command

kron(speye(10),A) than the size is 100x90. If matrix A would be square than a multiplication like this would certainly work, but alas it is not.

If anyone has a clue on how to do this I would love to hear it.

Regards,

Vincent


All times are GMT -4. The time now is 23:18.