3d polar plot for numerical solution of PDE system
1 view (last 30 days)
Show older comments
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 Comment
Bruno Pop-Stefanov
on 20 Jan 2014
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?
Answers (2)
Bruno Pop-Stefanov
on 20 Jan 2014
Edited: Bruno Pop-Stefanov
on 20 Jan 2014
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 .
0 Comments
Danila Zharenkov
on 21 Jan 2014
Edited: Danila Zharenkov
on 21 Jan 2014
1 Comment
Bruno Pop-Stefanov
on 21 Jan 2014
Edited: Bruno Pop-Stefanov
on 21 Jan 2014
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;
See Also
Categories
Find more on Polar Plots in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!