3d polar plot for numerical solution of PDE system
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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.
![](https://www.mathworks.com/matlabcentral/images/broken_image.png)
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 commentaire
Bruno Pop-Stefanov
le 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?
Réponses (2)
Bruno Pop-Stefanov
le 20 Jan 2014
Modifié(e) : Bruno Pop-Stefanov
le 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 commentaires
Danila Zharenkov
le 21 Jan 2014
Modifié(e) : Danila Zharenkov
le 21 Jan 2014
1 commentaire
Bruno Pop-Stefanov
le 21 Jan 2014
Modifié(e) : Bruno Pop-Stefanov
le 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;
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/177464/image.png)
Voir également
Catégories
En savoir plus sur Polar Plots dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!