MATLAB: Error with system of parabolic pdes

Question

Hi I am trying to solve the equations du/dt=-a1*u + f(u,v), dv/dt=-a2*v + g(u,v)

For simplicity, I now just take f and g as scalars

However I get the error message

Error using .* Matrix dimensions must agree.

Error in pdeasma (line 9) aod=a.*ar/12; % Off diagonal element

Error in assema (line 189) km1=pdeasma(it1,it2,it3,np,ar,x,y,sd,u,ux,uy,time,a(k,:));

Error in parabolic (line 84) [unused,MM]=assema(p,t,0,d,f,time);

Error in twoDpdex4 (line 26) u1=parabolic(u0,tlist,'squareb1',p,e,t,c,a,f,d);

The code is

c=zeros(2,2,2,2);

a=[-a1 0; 0 -a2];

d=[1 0; 0 1];

f=[1;2];

[p,e,t]=initmesh('squareg');

u0=zeros(size(p,2),2);

ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);

u0(ix,1)=0.1*ones(size(ix));

u0(ix,2)=0.2*ones(size(ix));

tlist=linspace(0,0.05,10);

u1=parabolic(u0,tlist,'squareb1',p,e,t,c,a,f,d);

Thank you in advance for your help and advice

Answer 1

Hi,

I'm afraid you've run into a couple of PDE Toolbox idiosyncrasies.

First, all coefficients must be column vectors (size n x 1 where various options are allowed for n). See this page for the c-coefficient, for example: http://www.mathworks.com/help/pde/ug/c.html For your case, if you want c to be all zero, a scalar, c=0, will work. A simple conversion of your d and a to column vectors fixes these.

Your definition of u0 as a number_of_points x N matrix makes definition easier. But PDE Toolbox requires this also to be a column vector.

The "squareb" boundary function example supplied with PDE Toolbox applies only to a single equation (scalar case). A different definition is required for your two-equation system. I supplied one that sets the values of both variables to zero on all edges. You can find more interesting examples here: http://www.mathworks.com/help/pde/ug/boundary-conditions-for-pde-systems.html

I've modified your code as shown below. Your code doesn't define a1 and a2 so I just took both of them to be one in my modified version.

Bill

----------------------------------------------------------------------------

function cmatrix_isaac(  )
%CMATRIX_ISSAC Summary of this function goes here
%   Detailed explanation goes here
c=0;
a1 = 1; a2 = 1;
a=[-a1 0 0 -a2]';
d=[1 0 0 1]';
f=[1;2];
[p,e,t]=initmesh('squareg');
u0=zeros(size(p,2),2);
ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);
u0(ix,1)=0.1*ones(size(ix));
u0(ix,2)=0.2*ones(size(ix));
tlist=linspace(0,0.05,10);
b = @boundaryFileZeroDirichlet;
u=parabolic(u0(:),tlist,b,p,e,t,c,a,f,d);
numTimes = size(u,2);
u2 = reshape(u, [], 2, numTimes);
for i=1:2
figure; pdeplot(p,e,t,'xydata', u2(:, i, end), 'contour', 'on');
end
end
function [ q, g, h, r ] = boundaryFileZeroDirichlet( p, e, u, time )
%BOUNDARYFILEZERODIRICHLET Solution is zero on all edges
%   Define a Dirichlet boundary condition making the solution equal zero
%   on all boundary edges.
N=2; % two-equation system
ne = size(e,2); % number of boundary edges
q = zeros(N^2, ne); % Neumann coefficient q is zero on all edges
g = zeros(N, ne); % Neumann coefficient g is zero on all edges
iN = eye(N); iN = iN(:);
h = iN*ones(1, 2*ne); % Dirichlet h coefficient is one at both ends of all edges
r = zeros(N,2*ne); % Dirichlet r coefficient is zero at both ends of all edges
end