How can I go about solving the following PDE?
2 次查看(过去 30 天)
显示 更早的评论
Hey, guys
Below is the equation I am trying to solve:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/151672/image.png)
I spent days trying to find appropriate solution adhering to the physical nature of problem. The desired result p(x,t) always give the initial value, and doesn't vary with distance although it should.
Can you give some hints or advice how to go forward? I need fresh advice.
Main file diffusivity_equation.m
*%
clear all
clc
%Parameters shared with the ODE routine
global n ncall phi k c cf ct mu h rw q beta xu xl re B ncase xw IC dx p0 pi
phi=0.2;
p0=3000;
k=125;
c=1.0e-08;
B=1.1;
cf=1.0e-06;
ct=c+cf;
mu=1.0;
h=100;
xw=3;
q=10000;
beta=-(q*B*mu/0.001127*k*h);
IC=-(q*mu)/(k*h*xw);
n=41;
xl=0;
xu=2000;
dx=(xu-xl)/(n-1);
for i=1:n
x(i)=(i-i)*dx;
end
%
%Independent variable for ODE integration
t0=0.0;
tf=12;
tout=linspace(t0,tf,n);
nout=n;
ncall=0;
%
p0=3000*ones(n,1);
%ODE integration
mf=1;
reltol=1.0e-04; abstol=1.0e-04;
options=odeset('RelTol',reltol,'AbsTol',abstol);
if(mf==1)%FDs
[t,p]=ode15s(@function_handle2,tout,p0,options);end
% Display output
for it=1:nout
fprintf('\n t = %4.2f\n',t(it));
for i=1:n
fprintf(' x = %4.1f p(x,t) = %8.5f\n',x(i),p(it,i));
end
end
fprintf('\n ncall = %5d\n',ncall);
%
ncall=ncall+1;*
function_hadle2.m
function pt=function_handle2(t,p)
*global n ncall phi k c cf ct mu h rw q beta xu xl re B px pxx ndss ncase
xw IC dx
%
%Problem parameters
xl=0.0;
xu=2000;
%
%left boundary condition
%
p(1)=0;
%
%first-order derivative px
%
px=dss004(0,xu,n,p)
%
%right boundary condition
px(n)=0;
%second-order dervative pxx
%
pxx=dss004(0,xu,n,px);
for i=1:n
pt(i)=beta*(cf*px(i)^2+pxx(i))
end
%PDE
pt=pt';
%
%Increment calls to function_handle2
ncall=ncall+1; *
dss004
% File: dss004.m
%
function [ux]=dss004(xl,xu,n,u)
ux( 2)=r4fdx*...
( -3.*u( 1) -10.*u( 2) +18.*u( 3) -6.*u( 4) +1.*u( 5));
for i=3:nm2
ux( i)=r4fdx*...
( +1.*u(i-2) -8.*u(i-1) +0.*u( i) +8.*u(i+1) -1.*u(i+2));
end
ux(n-1)=r4fdx*...
( -1.*u(n-4) +6.*u(n-3) -18.*u(n-2) +10.*u(n-1) +3.*u( n));
ux( n)=r4fdx*...
( 3.*u(n-4) -16.*u(n-3) +36.*u(n-2) -48.*u(n-1) +25.*u( n));
0 个评论
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!