Surface plot of PDE numeric solution

Hello, I'm solving the system of 2 PDE's, each function depends on 3 variables(radius, angle and time). I'm using explicit difference scheme. As a result I've got two 3D matrices (one for each function). To visualise the results, I want to plot surface plots for each function with fixed t(time). For example: In a moment t=0.5 the surface for u(:,:,0.5): X-axis will be the radius, Y-axis will be the angle and Z-axis will be the function value in in this point (x,y). Will be glad for any advices. Thanks.

Here is the code, if it's difficult to understand my English.

if true
  % code
   clc;
clear all;
%grid for r
n=100;
r_min=0.01;
r_max=1;
hr=(r_max-r_min)/(n-1)
r=r_min:hr:r_max
nr=max(size(r))
%grid for theta
m=100;
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
l=10000;
t_min=0;
t_max=1;
ht=(t_max-t_min)/(l-1)
time=t_min:ht:t_max
nt=max(size(time))
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
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
v(:,:,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;
            dudr=(u(i+1,j,t)-u(i-1,j,t))/(2*hr);
            d2udth2=(u(i,j+1,t)-2*u(i,j,t)+u(i,j-1,t))/hth^2;
            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,j,t))/hr;
            d2vdth2=(v(i,j+1,t)-2*v(i,j,t)+v(i,j-1,t))/hth^2;
            %centre nodes
            u(i,j,t+1)=u(i,j,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*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+(1/r(i)^2)*d2vdth2);%+a*v(i,j,t)*(1-c*u(i,j,t)-b*v(i,j,t)));

        end
        %derivatives for left boundary
        d2udr2= (u(i+1,1,t)-2*u(i,1,t)+u(i-1,1,t))/hr^2;
        dudr=(u(i+1,1,t)-u(i-1,1,t))/(2*hr);
        d2udth2=(u(i,2,t)-2*u(i,1,t)+u(i,nth,t))/hth^2;
        d2vdr2= (v(i+1,1,t)-2*v(i,1,t)+v(i-1,1,t))/hr^2;
        dvdr=(v(i+1,1,t)-v(i-1,1,t))/(2*hr);
        d2vdth2=(v(i,2,t)-2*v(i,1,t)+v(i,nth,t))/hth^2;        
        uleft=u(i,1,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*d2udth2);%+u(i,1,t)*(1-u(i,1,t)-c*v(i,1,t)));
        vleft=v(i,1,t)+ht*mu*(d2vdr2+(1/r(i))*dvdr+(1/r(i)^2)*d2vdth2);%+a*v(i,1,t)*(1-c*u(i,1,t)-b*v(i,1,t)));        
        %derivatives for right boundary
        d2udr2= (u(i+1,nth,t)-2*u(i,nth,t)+u(i-1,nth,t))/hr^2;
        dudr=(u(i+1,nth,t)-u(i-1,nth,t))/(2*hr);
        d2udth2=(u(i,1,t)-2*u(i,nth,t)+u(i,nth-1,t))/hth^2;
        d2vdr2= (v(i+1,nth,t)-2*v(i,nth,t)+v(i-1,nth,t))/hr^2;
        dvdr=(v(i+1,nth,t)-v(i-1,nth,t))/(2*hr);
        d2vdth2=(v(i,1,t)-2*v(i,nth,t)+v(i,nth-1,t))/hth^2;        
        uright=u(i,nth,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*d2udth2);%+u(i,nth,t)*(1-u(i,nth,t)-c*v(i,nth,t)));
        vright=v(i,nth,t)+ht*mu*(d2vdr2+(1/r(i))*dvdr+(1/r(i)^2)*d2vdth2);%+a*v(i,nth,t)*(1-c*u(i,nth,t)-b*v(i,nth,t))); 
        %filling boundaries for theta
        u(i,1,t+1)=uleft;
        v(i,1,t+1)=vleft;
        u(i,nth,t+1)=uright;
        v(i,nth,t+1)=vright;
    end
    %boundaries for r
    for j=1:nth
        u(1,j,t+1)=u(2,j,t+1);
        u(nr,j,t+1)=u(nr-1,j,t+1);
        v(1,j,t+1)=v(2,j,t+1);
        v(nr,j,t+1)=v(nr-1,j,t+1);
    end
end
  end

surf(r,theta,?)

[T K]

Could you please send me a picture of the second degree system of equations with boundary conditions?

Answer 1

[Youssef Khmou]

Danila, I tried to examine you solution, the first thing is that r0,r1,r2 and r3 are not provided , plus a variable mu is not defined, i believe mu is an internal variable name, so i changed it to Mu to make difference, here is the changed version of your solution. Concerning the last question yo can simply use surf(r,theta,u(:,:,n)) for the first instant n=1, n=2,..... n=end :

 clc;clear all;
n=100;
r_min=0.01;
r0=r_min;
r_max=1;
r3=r_max;
r1=r0+0.2;
r2=r1+0.2;
hr=(r_max-r_min)/(n-1);
r=r_min:hr:r_max;
nr=max(size(r));
%grid for theta
m=100;
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
l=100;
t_min=0;
t_max=1;
ht=(t_max-t_min)/(l-1);
time=t_min:ht:t_max;
nt=max(size(time));
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 || r(i)>r3)
          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);
      end
      u(i,j,1)=u0;
  end
end
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
v(:,:,1);
Mu=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);
            dudr=(u(i+1,j,t)-u(i-1,j,t))/(2*hr);
            d2udth2=(u(i,j+1,t)-2*u(i,j,t)+u(i,j-1,t))/(hth^2);
            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,j,t))/hr;
            d2vdth2=(v(i,j+1,t)-2*v(i,j,t)+v(i,j-1,t))/(hth^2);
            %centre nodes
            u(i,j,t+1)=u(i,j,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*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+(1/r(i)^2)*d2vdth2);%+a*v(i,j,t)*(1-c*u(i,j,t)-b*v(i,j,t)));
        end
        %derivatives for left boundary
        d2udr2= (u(i+1,1,t)-2*u(i,1,t)+u(i-1,1,t))/hr^2;
        dudr=(u(i+1,1,t)-u(i-1,1,t))/(2*hr);
        d2udth2=(u(i,2,t)-2*u(i,1,t)+u(i,nth,t))/hth^2;
        d2vdr2= (v(i+1,1,t)-2*v(i,1,t)+v(i-1,1,t))/hr^2;
        dvdr=(v(i+1,1,t)-v(i-1,1,t))/(2*hr);
        d2vdth2=(v(i,2,t)-2*v(i,1,t)+v(i,nth,t))/hth^2;        
        uleft=u(i,1,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*d2udth2);%+u(i,1,t)*(1-u(i,1,t)-c*v(i,1,t)));
        vleft=v(i,1,t)+ht*Mu*(d2vdr2+(1/r(i))*dvdr+(1/r(i)^2)*d2vdth2);%+a*v(i,1,t)*(1-c*u(i,1,t)-b*v(i,1,t)));        
        %derivatives for right boundary
        d2udr2= (u(i+1,nth,t)-2*u(i,nth,t)+u(i-1,nth,t))/hr^2;
        dudr=(u(i+1,nth,t)-u(i-1,nth,t))/(2*hr);
        d2udth2=(u(i,1,t)-2*u(i,nth,t)+u(i,nth-1,t))/hth^2;
        d2vdr2= (v(i+1,nth,t)-2*v(i,nth,t)+v(i-1,nth,t))/hr^2;
        dvdr=(v(i+1,nth,t)-v(i-1,nth,t))/(2*hr);
        d2vdth2=(v(i,1,t)-2*v(i,nth,t)+v(i,nth-1,t))/hth^2;        
        uright=u(i,nth,t)+ht*(d2udr2+(1/r(i))*dudr+(1/r(i)^2)*d2udth2);%+u(i,nth,t)*(1-u(i,nth,t)-c*v(i,nth,t)));
        vright=v(i,nth,t)+ht*Mu*(d2vdr2+(1/r(i))*dvdr+(1/r(i)^2)*d2vdth2);%+a*v(i,nth,t)*(1-c*u(i,nth,t)-b*v(i,nth,t))); 
        %filling boundaries for theta
        u(i,1,t+1)=uleft;
        v(i,1,t+1)=vleft;
        u(i,nth,t+1)=uright;
        v(i,nth,t+1)=vright;
    end
    %boundaries for r
    for j=1:nth
        u(1,j,t+1)=u(2,j,t+1);
        u(nr,j,t+1)=u(nr-1,j,t+1);
        v(1,j,t+1)=v(2,j,t+1);
        v(nr,j,t+1)=v(nr-1,j,t+1);
    end
end

 figure, surf(r,theta,u(:,:,3)), shading interp, title(' U_{3}');
 figure, surf(r,theta,v(:,:,3)), shading interp, title(' V_{3}');

Answer 2

[Danila Zharenkov]

Thank you. The constants are defined, it's okay, I just remove them to reduce code on forum. This works, thanks again.