How to stabilize discretized ODEs by method of lines?
2 次查看(过去 30 天)
显示 更早的评论
Hello, guys
I am about to solve linear diffusivity equation representing pressure distribution in reservoir for transient conditions. Originally, the pressure should gradually i.e. exponentially rise as it goes further from the center of the well. However, matlab shows pressure declining tendency as it goes from the well. How can I stabilize the solution so that I can obtain proper solution? The correct solution should be pressure gradually increasing from x=0 to x=1000.
The codes are attached below:
*% lineardiffusion1.m*
global ncall D x0 IP
D=3000;
x0=1000;
n=51;
dx=x0/(n-1);
IP=3500;
u0=IP*ones(n,1);
for i=1:n
x(i)=(i-1)*dx+50;
end
tf=50
tout=[0.0:2:tf]';
nout=length(tout);
ncall=0;
reltol=1.0e-06; abstol=1.0e-06;
options=odeset('RelTol',reltol,'AbsTol',abstol);
[t,u]=ode15s(@lineardiffusion,tout,u0,options)
%
% Display a heading and selected output
fprintf('\n nx = %2d x0 = %4.2f std = %4.2f\n',nx,x0,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:nx
fprintf(' x = %4.1f u(x,t) = %8.5f\n',x(i),u(it,i));
end
end
fprintf('\n ncall = %5d\n',ncall);
%
% Parametric plot
figure(1);
plot(x,u); axis([0 x0 0 10000]);
title('u(x,t), t = 0, 0.25, ..., 2.5'); xlabel('x'); ylabel('u(x,t)')
surf(x, t ,u);
%contour(r, tt, u2);
axis([0 x0 0 tf 0 10000]);
xlabel('Linear 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,x,u(1,:),'k.')
time(1)=0;
hold on
plot(axes2,x,u(10,:),'m.')
time(2)=t(10);
plot(x,u(15,:),'b.')
time(3)=t(15);
plot(x,u(20,:),'g.')
time(4)=t(20);
plot(x,u(30,:),'c.')
time(4)=t(30);
plot(x,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 'doc.pdf'
*%lineardiffusion*
function ut=lineardiffusion(t,u)
global ncall D x0 IP xl xu nx
nx=51;
for i=1:nx
ux=dss008(0.0,x0,nx,u);
u(nx)=IP;
uxx=dss008(0.0,x0,nx,ux);
ux(1)=-100;
ut(i)=D*(uxx(i));
end
ut=ut';
ut(nx)=0;
ncall=ncall+1;
% File: dss008.m
%
function [ux]=dss008(xl,xu,n,u)
% Compute the spatial increment
dx=(xu-xl)/(n-1);
r8fdx=1./(40320.*dx);
nm4=n-4;
%
% Grid point 1
ux( 1)=r8fdx*...
(-109584. *u( 1)...
+322560. *u( 2)...
-564480. *u( 3)...
+752640. *u( 4)...
-705600. *u( 5)...
+451584. *u( 6)...
-188160. *u( 7)...
+46080. *u( 8)...
-5040. *u( 9));
%
% Grid point 2
ux( 2)=r8fdx*...
(-5040. *u( 1)...
-64224. *u( 2)...
+141120. *u( 3)...
-141120. *u( 4)...
+117600. *u( 5)...
-70560. *u( 6)...
+28224. *u( 7)...
-6720. *u( 8)...
+720. *u( 9));
%
% Grid point 3
ux( 3)=r8fdx*...
(+720. *u( 1)...
-11520. *u( 2)...
-38304. *u( 3)...
+80640. *u( 4)...
-50400. *u( 5)...
+26880. *u( 6)...
-10080. *u( 7)...
+2304. *u( 8)...
-240. *u( 9));
%
% Grid point 4
ux( 4)=r8fdx*...
(-240. *u( 1)...
+2880. *u( 2)...
-20160. *u( 3)...
-18144. *u( 4)...
+50400. *u( 5)...
-20160. *u( 6)...
+6720. *u( 7)...
-1440. *u( 8)...
+144. *u( 9));
%
% Grid point i
for i=5:nm4
ux( i)=r8fdx*...
(+144. *u(i-4)...
-1536. *u(i-3)...
+8064. *u(i-2)...
-32256. *u(i-1)...
+0. *u(i )...
+32256. *u(i+1)...
-8064. *u(i+2)...
+1536. *u(i+3)...
-144. *u(i+4));
end
%
% Grid point n-3
ux(n-3)=r8fdx*...
(-144. *u(n-8)...
+1440. *u(n-7)...
-6720. *u(n-6)...
+20160. *u(n-5)...
-50400. *u(n-4)...
+18144. *u(n-3)...
+20160. *u(n-2)...
-2880. *u(n-1)...
+240. *u(n ));
%
% Grid point n-2
ux(n-2)=r8fdx*...
(+240. *u(n-8)...
-2304. *u(n-7)...
+10080. *u(n-6)...
-26880. *u(n-5)...
+50400. *u(n-4)...
-80640. *u(n-3)...
+38304. *u(n-2)...
+11520. *u(n-1)...
-720. *u(n ));
%
% Grid point n-1
ux(n-1)=r8fdx*...
(-720. *u(n-8)...
+6720. *u(n-7)...
-28224. *u(n-6)...
+70560. *u(n-5)...
-117600. *u(n-4)...
+141120. *u(n-3)...
-141120. *u(n-2)...
+64224. *u(n-1)...
+5040. *u(n ));
%
% Grid point n
ux(n )=r8fdx*...
(+5040. *u(n-8)...
-46080. *u(n-7)...
+188160. *u(n-6)...
-451584. *u(n-5)...
+705600. *u(n-4)...
-752640. *u(n-3)...
+564480. *u(n-2)...
-322560. *u(n-1)...
+109584. *u(n ));
0 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Boundary Conditions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!