Error with system of parabolic pdes
5 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
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
More Answers (2)
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
0 Comments
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
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!