Error with system of parabolic pdes

5 views (last 30 days)
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

Accepted Answer

Bill Greene
Bill Greene on 15 May 2013
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
  1 Comment
Isaac
Isaac on 16 May 2013
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

Sign in to comment.

More Answers (2)

Bill Greene
Bill Greene on 16 May 2013
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

Bill Greene
Bill Greene on 16 May 2013
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
  1 Comment
Isaac
Isaac on 16 May 2013
A tremendous thank you to you, Sir. I will upgrade my Matlab and hopefully have this all sorted out.
Thanks again, Isaac

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!