% File: pde_1_main_3D.m % From Chapter 14 - 'Diffusion Equation in Spherical Coordinates', % of the book: % William E Schiesser and Graham W Griffiths (2009). % A Compendium of Partial Differential Equation Models % - Method of Lines Analysis with Matlab, % Cambridge University Press (ISBN-13: 9780521519861). % % Produces the 3D plots % % Clear previous files clear all clc close all % % Global area global r r0 nr D tau ncall IP %f % % Model parameters D=1000; r0=1000; std=1.0; tau=1.0; IP=2500; %initial pressure % % Radial grid and inhomogeneous term nr=51; dr=r0/(nr-1); for i=1:nr r(i)=(i-1)*dr+200; %f(i)=exp(-r(i)^2/std^2); %f(i)=-0.000005; %set constant mass extraction rate; end %f(1:nr)=0.0; %f(1)=-0.25 % % Independent variable for ODE integration tf= 50; tau=tf; tout=[0.0:2:tf]'; nout=length(tout); ncall=0; % % Initial condition % u0=zeros(nr,1); u0=IP*ones(nr,1); % % ODE integration reltol=1.0e-06; abstol=1.0e-06; options=odeset('RelTol',reltol,'AbsTol',abstol); [t,u]=ode15s(@pde_1_oil,tout,u0,options); % % Display a heading and selected output fprintf('\n nr = %2d r0 = %4.2f std = %4.2f\n',nr,r0,std); disp('at one') % pause for it=1:nout it fprintf('\n t = %4.2f\n',t(it)); disp('at two') %pause for i=1:2:nr fprintf(' r = %4.1f u(r,t) = %8.5f\n',r(i),u(it,i)); end end fprintf('\n ncall = %5d\n',ncall); % % Parametric plot figure(1); plot(r,u); axis([0 r0 0 3000]); title('u(r,t), t = 0, 0.25, ..., 2.5'); xlabel('r'); ylabel('u(r,t)') % % Extend integration for closer approach to steady state fprintf('\n Extended solution\n'); u0=u(nout,:); t0=t(nout); tf=t0+5.0; tout=[t0:0.25:tf]'; nout=length(tout); [t1,u1]=ode15s(@pde_1_oil,tout,u0,options); for it=1:nout fprintf('\n t = %4.1f\n',t1(it)); for i=1:2:nr fprintf(' r = %4.1f u(r,t) = %8.5f\n',r(i),u1(it,i)); end end fprintf('\n ncall = %5d\n',ncall); tt=[t;t1]; u2=[u;u1]; surf(r, tt ,u2); %contour(r, tt, u2); axis([0 r0 0 tf 2000 2700]); xlabel('radial distance'); ylabel('time'); zlabel('pressure'); title([{'Surface plot at time t=',tf}])%, ... % { 'at time t=7.5'}]); rotate3d on disp('check plot and press return to progress') pause hold 'off' % Create figure figure2 = figure('PaperSize',[20.98 29.68]); % Create axes axes2 = axes('Parent',figure2); box('on'); hold('all'); si=size(u,1); plot(axes2,r,u(1,:),'k.') time(1)=0; hold on plot(axes2,r,u(1,:),'m.') time(2)=t(1); plot(r,u(5,:),'b.') time(3)=t(5); plot(r,u(10,:),'g.') time(4)=t(10); plot(r,u(15,:),'c.') time(4)=t(15); plot(r,u(si,:),'r.') time(5)=t(si); % Create legend legend2=legend(axes2,'t=0','t=0.225','t=0.625','t=1.8725','t=2.5'); % this works set(legend2,'Position',[0.5354 0.2781 0.1946 0.2333]); print -dpdf 'forDani.pdf'