3d polar plot for numerical solution of PDE system

1 次查看（过去 30 天）
Danila Zharenkov
Danila Zharenkov 2014-1-20
编辑： Bruno Pop-Stefanov 2014-1-21
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"
  1 个评论
Bruno Pop-Stefanov
Bruno Pop-Stefanov 2014-1-20
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?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Bruno Pop-Stefanov
Bruno Pop-Stefanov 2014-1-20
编辑：Bruno Pop-Stefanov 2014-1-20
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 .
Danila Zharenkov
Danila Zharenkov 2014-1-21
编辑：Danila Zharenkov 2014-1-21
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
  1 个评论
Bruno Pop-Stefanov
Bruno Pop-Stefanov 2014-1-21
编辑：Bruno Pop-Stefanov 2014-1-21
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;

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 PDE Solvers 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by