How can I enter a source term which is a vector, to solve -div(cgradu)=F(F a vector)?

Hello, I want to solve the pde (-div(cgrad(u))= F) with F given as a column vector of length N. I used the assempde to solve directly the pde with that F but it always answers me

??? Error using ==> plus Matrix dimensions must agree.

Error in ==> assempde at 245 KK=K+M+Q;

Error in ==> myprogram at 119 u=assempde(b,p,e,t,c,a,F);

Any help, please? Thank you in advance

[Walter Roberson]

At the MATLAB command line, command

 dbstop if error

and then run the program. When it stops, what does size(K), size(M), size(Q) indicate? Then, give the command

dbup

and show size(b), size(p), size(e), size(t), size(c), size(a), size(F)

Answer 1

[Bill Greene]

I assume N is > 1? How did you define the first argument to assempde, b? I'm guessing that the problem at line 245 of assempde is caused by K and M having a different size from Q where Q is calculated based on your boundary condition definition, b.

Bill

[Batoul]

yes, N is strictly greater than 1 (it is equal to the number of nodes) and for b, I defined it using pdebound command function.

[Bill Greene]

OK, well I suspect there is a problem in your pdebound function. There is no way to say more with only the information you have provided.

Bill

[Batoul]

Yes but when i give F as a constant(i.e. N=1), it works... and for the pdebound i used this

   function [qmatrix,gmatrix,hmatrix,rmatrix] = pdebound(p,e,~,~)
  ne = size(e,2); % number of edges

qmatrix = zeros(1,ne);
gmatrix = qmatrix;
hmatrix = zeros(1,2*ne);
rmatrix = hmatrix;

for k = 1:ne
            hmatrix(k) = 1;
            hmatrix(k+ne) = 1;
            rmatrix(k) =0;
            rmatrix(k+ne) =0;

end

since i have u=0 on the boundary

So, i can't find the mistake :(

Answer 2

[Bill Greene]

Yes, your pdebound function is incorrect for a system of PDE (N>1). Specifically, as this documentation page shows,

http://www.mathworks.com/help/pde/ug/boundary-conditions-for-pde-systems.html

the returned matrices must have these dimensions (for this example, N=3):

N = 3; % Set N = the number of equations
ne = size(e,2); % number of edges
qmatrix = zeros(N^2,ne);
gmatrix = zeros(N,ne);
hmatrix = zeros(N^2,2*ne);
rmatrix = zeros(N,2*ne);

I recommend taking a close look at the www page I list above so you can see exactly how those matrices must be defined.

Bill