Error with system of parabolic pdes

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

[Bill Greene]

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

[Isaac]

Hi Bill,

Thank you so much for your last response. It really helped me out. However, I wonder if it is possible to pick your brains a little more on this topic. I now wish to call f in a functional form ie

 function f = fcoeffunction(p,t,u,time)

    N = 2; % Number of equations
    % Triangle point indices
    it1=t(1,:);
    it2=t(2,:);
    it3=t(3,:);

    % Find centroids of triangles
    xpts=(p(1,it1)+p(1,it2)+p(1,it3))/3;
    ypts=(p(2,it1)+p(2,it2)+p(2,it3))/3;

    [ux,uy] = pdegrad(p,t,u); % Approximate derivatives
    uintrp = pdeintrp(p,t,u); % Interpolated values at centroids
    uintrp = reshape(uintrp,[],N); % matrix with N column
    uintrp = uintrp'; % change to row vectors

    nt = size(t,2); % Number of columns
    f = zeros(N,nt); % Allocate f

     Now the particular functional form of f
    f(1,:) = xpts - ypts + uintrp(1,:);
    f(2,:) = 1 + tanh(ux(1,:)) + tanh(uy(3,:));
end

The rest of the code that you kindly provided is still the same. However, when I try to call it using either

 1 f = @fcoeffunction;
 u=parabolic(u0(:),tlist,b,p,e,t,c,a,f,d);

or

 2 u=parabolic(u0(:),tlist,b,p,e,t,c,a,'fcoeffunction',d);

I get the error

 Error using assema (line 223)
Wrong number of rows of a.

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

 Error in twopde (line 21)
    u=parabolic(u0(:),tlist,b,p,e,t,c,a,'fcoeffunction',d);

Could you please give me some advice on the best way to go about fixing this issue?

Thanks again, Isaac

Answer 2

[Bill Greene]

Hi,

You've run into a PDE Toolbox bug that has been fixed for several releases. Is there any possibility that you can upgrade to the latest version of MATLAB (R2013a)?

If you are able to upgrade, a couple of small changes in your fcoeffunction will get this working. First, these two lines should be removed:

uintrp = reshape(uintrp,[],N); % matrix with N column
uintrp = uintrp'; % change to row vectors

Then, this line

f(2,:) = 1 + tanh(ux(1,:)) + tanh(uy(3,:));

is referring to uy(3,:) which doesn't exist for a 2-equation system.

Regards,

Bill

Answer 3

[Bill Greene]

I should have also mentioned that the functionality to have coefficients that are functions of the solution (or gradients) in the parabolic function was a new capability in PDE Toolbox R2012b. So, if you are relying on that capability, you will need to upgrade to at least R2012b.

Bill

[Isaac]

A tremendous thank you to you, Sir. I will upgrade my Matlab and hopefully have this all sorted out.

Thanks again, Isaac