3d polar plot for numerical solution of PDE system

Hello, I;m solving the system of PDE in polar coordinates on a disk. Results must be present as a 3d polar plot, something like on the attached picture.

I've got the vector of Theta, vector of radiuses and a matrix of values for respected Theta and R.

Here is the code:

      clc;
  clear all;
  filename = 'plant.gif';
  %grid for r
  n=30;
  r_min=0.01;
  r_max=1;
  hr=(r_max-r_min)/(n-1);
  r(1)=r_min-hr;
  r(2:n+1)=r_min:hr:r_max;
  r(n+2)=r(n+1)+hr;
  nr=max(size(r));

%grid for theta
m=30;
th_min=0;
th_max=2*pi;
hth=(th_max-th_min)/(m-1);
theta=th_min:hth:th_max;
nth=max(size(theta));

%grid for t
t_min=0;
t_max=1;
l=(2*(n+2)^2)*t_max+1;
ht=(t_max-t_min)/(l-1);
time=t_min:ht:t_max;
nt=max(size(time))

%constants
a=0.99;
b=1;
c=100;
mu=1;
r0=0.2;
r1=0.4;
r2=0.6;
r3=0.8;

u = zeros(nr,nth,nt);
v = zeros(nr,nth,nt);

%Inititalization
u0=0;
for i=1:nr
    for j=1:nth
        if (r(i)<r0)
            u0=0;
        elseif ((r(i)>=r0) &&(r(i)<r1))
            u0=(r(i)-r0)/(r1-r0);
        elseif ((r(i)>=r1)&&(r(i)<r2))
            u0=1;
        elseif ((r(i)>=r2)&&(r(i)<r3))
            u0=(r(i)-r3)/(r2-r3);
        else
            u0=0;
        end
        u(i,j,1)=u0;
    end
end
[bf,af]=butter(2,0.5);
u(:,:,1)=filtfilt(bf,af,u(:,:,1));

v0=0;
for i=1:nr
    for j=1:nth
        if(r(i)<r1)
            v0=1;
        elseif ((r(i)>=r1)&&(r(i)<r3))
            v0=(r(i)-r3)/(r1-r3);
        else
            v0=0;
        end
        v(i,j,1)=v0;
    end
end
[bf,af]=butter(2,0.5);
v(:,:,1)=filtfilt(bf,af,v(:,:,1));

S.fh = figure('Units','pixels', ...
                     'Name','EMGGUI', ...
                     'NumberTitle','off', ...
                     'MenuBar','none',...
                     'Toolbar','none',...
                     'Position',[200 100 950 500], ... %left bottom width height
                     'Resize','on', ...
                     'Visible','off');

S.ax(1) = axes('units','pixels',...
            'position',[50 50 400 400],...
            'drawmode','fast');

S.ax(2) = axes('units','pixels',...
            'position',[500 50 400 400],...
            'drawmode','fast');

axes(S.ax(1))
surf(theta,r,u(:,:,1))
        set(gca,'xlim',[0 2*pi])
        set(gca,'ylim',[0 1])
axes(S.ax(2))
surf(theta,r,v(:,:,1))
        set(gca,'xlim',[0 2*pi])
        set(gca,'ylim',[0 1])

for t=1:nt-1

      for i=2:nr-1
          for j=2:nth-1
              %derivatives

              d2udr2= (u(i+1,j,t)-2*u(i,j,t)+u(i-1,j,t))/hr^2;
              d2udth2= (u(i,j+1,t)-2*u(i,j,t)+u(i,j-1,t))/hr^2;
              dudr=(u(i+1,j,t)-u(i-1,j,t))/(2*hr);
              d2vdr2= (v(i+1,j,t)-2*v(i,j,t)+v(i-1,j,t))/hr^2;
              dvdr=(v(i+1,j,t)-v(i-1,j,t))/(2*hr);
              d2vdth2= (v(i,j+1,t)-2*v(i,j,t)+v(i,j-1,t))/hr^2;
              %centre nodes
              u(i,j,t+1)=u(i,j,t)+ht*(d2udr2+(1/r(i))*dudr+d2udth2+u(i,j,t)*(1-u(i,j,t)-c*v(i,j,t)));
              v(i,j,t+1)=v(i,j,t)+ht*mu*(d2vdr2+(1/r(i))*dvdr+d2vdth2+a*v(i,j,t)*(1-c*u(i,j,t)-b*v(i,j,t)));

          end
      end
      %boundaries for r
      for j=1:nth
          u(1,j,t+1)=u(3,j,t+1);
          u(nr,j,t+1)=u(nr-2,j,t+1);
          v(1,j,t+1)=v(3,j,t+1);
          v(nr,j,t+1)=v(nr-2,j,t+1);
      end
          u(2:end,1,t+1)=u(2:end,2,t+1);
          v(2:end,1,t+1)=v(2:end,2,t+1);
          u(2:end,nth,t+1)=u(2:end,nth-1,t+1);
          v(2:end,nth,t+1)=v(2:end,nth-1,t+1);
          %corners
          u(1,1,t+1)=u(1,2,t+1);
          u(1,nth,t+1)=u(1,nth-1,t+1);
          u(nr,1,t+1)=u(nr-1,1,t+1);
          u(nr,nth,t+1)=u(nr-1,nth,t+1);
          v(1,1,t+1)=v(1,2,t+1);
          v(1,nth,t+1)=v(1,nth-1,t+1);
          v(nr,1,t+1)=v(nr-1,1,t+1);
          v(nr,nth,t+1)=v(nr-1,nth,t+1);

axes(S.ax(1))
surf(theta,r,u(:,:,t))
xlabel('theta');
ylabel('r');
title('u');
        set(gca,'xlim',[0 2*pi])
        set(gca,'ylim',[0 1])
if max(u(:,1,t))>0.1
    set(gca,'zlim',[0 1])
elseif max(u(:,1,t))>0.01
    set(gca,'zlim',[0 0.1])
elseif max(u(:,1,t))>0.001
    set(gca,'zlim',[0 0.01])
end
axes(S.ax(2)) %#ok<LAXES>

surf(theta,r,v(:,:,t));
xlabel('theta');
ylabel('r');
title('v')
        set(gca,'xlim',[0 2*pi])
        set(gca,'ylim',[0 1])
if v(1,1,t)>0.1
    set(gca,'zlim',[0 1])
elseif v(1,1,t)>0.01
    set(gca,'zlim',[0 0.1])
elseif v(1,1,t)>0.001
    set(gca,'zlim',[0 0.01])
end
frame = getframe(S.fh);
      im = frame2im(frame);
      [imind,cm] = rgb2ind(im,256);
      if t == 1;
          imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
      else
          imwrite(imind,cm,filename,'gif','WriteMode','append');
      end
end

      % code
    end

I've tried to use polar2cart function and plot the surface, but matlab throws an error "Error using .* Matrix dimensions must agree"

[Bruno Pop-Stefanov]

I am trying to understand your code to plot what you would like to plot. If I am correct, you would like to plot u and v in 3D against theta and r?

Answer 1

[Bruno Pop-Stefanov]

When using the pol2cart function, theta, r and z must be of the same size. That is, your matrix of values for the corresponding theta and r should be a vector of the same size, and not a matrix. This is where your error is coming from; however, pol2cart lets you transform only cylindrical data, like a helix for example. This won't work for a surface plot.

To obtain a real 3D plot from polar coordinates like in the picture you are showing, you can use polar3d from the MATLAB Central File Exchange: http://www.mathworks.com/matlabcentral/fileexchange/7656-3d-polar-plot .

Answer 2

[Danila Zharenkov]

Thanks! But i've got a problem with this function

I'm trying to plot my figures, but it looks very strange. I don't now what to write in set() for polar plot. Can you give me some advice? Thanks.

Here's my code

clc;
clear all;

filename = 'plant.gif';

%grid for r
n=5;
r_min=0.01;
r_max=1;
hr=(r_max-r_min)/(n-1);
r(1)=r_min-hr;
r(2:n+1)=r_min:hr:r_max;
r(n+2)=r(n+1)+hr;
nr=max(size(r));

%grid for theta
m=5;
th_min=0;
th_max=2*pi;
hth=(th_max-th_min)/(m-1);
theta(2:m+1)=th_min:hth:th_max;
theta(1)=-hth;
nth=max(size(theta));

%grid for t
t_min=0;
t_max=1;
l=(2*(m+2)^2)*(2*(n+2)^2)*t_max+1;
ht=(t_max-t_min)/(l-1);
time=t_min:ht:t_max;
nt=max(size(time));

%constants
a=0.99;
b=1;
c=100;
mu=1;
r0=0.2;
r1=0.4;
r2=0.6;
r3=0.8;

u = zeros(nr,nth,nt);
v = zeros(nr,nth,nt);

for i=1:nr
    for j=1:nth
        if ((r(i)<=r0)&&(0<=theta(j))&&(theta(j)<=pi/6))
           u0=0; 
        elseif ((r(i)<=r0)&&(pi/6<theta(j))&&(theta(j)<2*pi))
           u0=1;
        elseif ((r(i)>r0)&&(0<=theta(j))&&(theta(j)<=2*pi))
           u0=0; 
        end
        if j==1 &&(r(i)<=r0)
            u0=1;
        end
        if ((r(i)<=r0)&&(pi/6<theta(j))&&(theta(j)==2*pi))
            u0=0;
        end
        u(i,j,1)=u0;
    end
end
% [bf,af]=butter(2,0.9);
% u(:,:,1)=filtfilt(bf,af,u(:,:,1));
v0=0;
for i=1:nr
    for j=1:nth
        if ((r(i)<=r0)&&(0<=theta(j))&&(theta(j)<=pi/6))
           v0=1; 
        elseif ((r(i)<=r0)&&(pi/6<theta(j))&&(theta(j)<2*pi))
           v0=0;
        elseif ((r(i)>r0)&&(0<=theta(j))&&(theta(j)<=2*pi))
           v0=0; 
        end
        if j==1
            v0=0;
        end
        if ((r(i)<=r0)&&(pi/6<theta(j))&&(theta(j)==2*pi))
            v0=1;
        end
        v(i,j,1)=v0;
    end
end
%  [bf,af]=butter(2,0.9);
%  v(:,:,1)=filtfilt(bf,af,v(:,:,1));

S.fh = figure('Units','pixels', ...
                     'Name','Virus Dynamcis', ...
                     'NumberTitle','off', ...
                     'MenuBar','none',...
                     'Toolbar','none',...
                     'Position',[200 100 950 500], ... %left bottom width height
                     'Resize','on', ...
                     'Visible','off');

S.ax(1) = axes('units','pixels',...
            'position',[50 50 400 400],...
            'drawmode','fast');

S.ax(2) = axes('units','pixels',...
            'position',[500 50 400 400],...
            'drawmode','fast');

axes(S.ax(1))
%surf(theta,r,u(:,:,1))
polar3d(u(:,:,1),0,2*pi,0,1,1,'surf');

axes(S.ax(2))
%surf(theta,r,v(:,:,1))
polar3d(v(:,:,1),0,2*pi,0,1,1,'surf');

for t=1:nt-1

      for i=2:nr-1
          for j=2:nth-1
              %derivatives

              d2udr2= (u(i+1,j,t)-2*u(i,j,t)+u(i-1,j,t))/hr^2;
              d2udth2= (u(i,j+1,t)-2*u(i,j,t)+u(i,j-1,t))/hr^2;
              dudr=(u(i+1,j,t)-u(i-1,j,t))/(2*hr);
              d2vdr2= (v(i+1,j,t)-2*v(i,j,t)+v(i-1,j,t))/hr^2;
              dvdr=(v(i+1,j,t)-v(i-1,j,t))/(2*hr);
              d2vdth2= (v(i,j+1,t)-2*v(i,j,t)+v(i,j-1,t))/hr^2;
              %centre nodes
              u(i,j,t+1)=u(i,j,t)+ht*(d2udr2+(1/r(i))*dudr+d2udth2+u(i,j,t)*(1-u(i,j,t)-c*v(i,j,t)));
              v(i,j,t+1)=v(i,j,t)+ht*mu*(d2vdr2+(1/r(i))*dvdr+d2vdth2+a*v(i,j,t)*(1-c*u(i,j,t)-b*v(i,j,t)));

          end

      end
      %boundaries for r
      for j=1:nth
          u(1,j,t+1)=u(3,j,t+1);
          u(nr,j,t+1)=u(nr-2,j,t+1);
          v(1,j,t+1)=v(3,j,t+1);
          v(nr,j,t+1)=v(nr-2,j,t+1);
      end
          u(:,1,t+1)=u(:,nth-1,t+1);
          v(:,1,t+1)=v(:,nth-1,t+1);
          u(:,nth,t+1)=u(:,2,t+1);
          v(:,nth,t+1)=v(:,2,t+1);
  %         %corners
          u(1,1,t+1)=u(1,2,t+1);
          u(1,nth,t+1)=u(1,nth-1,t+1);
          u(nr,1,t+1)=u(nr-1,1,t+1);
          u(nr,nth,t+1)=u(nr-1,nth,t+1);
          v(1,1,t+1)=v(1,2,t+1);
          v(1,nth,t+1)=v(1,nth-1,t+1);
          v(nr,1,t+1)=v(nr-1,1,t+1);
          v(nr,nth,t+1)=v(nr-1,nth,t+1);

axes(S.ax(1))
%surf(theta,r,u(:,:,t))
polar3d(u(:,:,t),0,2*pi,0,1,1,'surf');
title('u');

end
axes(S.ax(2))
polar3d(v(:,:,t),0,2*pi,0,1,1,'surf');
%surf(theta,r,v(:,:,t))
title('v')

frame = getframe(S.fh);
      im = frame2im(frame);
      [imind,cm] = rgb2ind(im,256);
      if t == 1;
          imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
      else
          imwrite(imind,cm,filename,'gif','WriteMode','append');
      end

end

[Bruno Pop-Stefanov]

Your input matrices u and v are very small: 7x6 at each time step. Try to discretize your space in theta and r more and you won't see disgracious planes in your plot. You might get an out-of-memory error because your discretization in time is so fine. Discretize t with a fixed step size so that it does not depend on m or n.

For example, try with:

n = 10;
m = 10;
ht = 0.001;