solving 1D parabolic-elliptic equations with one initial condition

Hello, I am asking for help for solving a system of partial differential equations on 1D. I have two equaitons, one is parabolic and another one is elliptic. But I have only one intial value for parabolic equation. I am trying to use pdepe solver. However, pdepe solver requires that the number of initial conditions must equal the number of equations. Is there any Matlab solver can slove such problem?

Answer 1

[Torsten]

Supply an arbitrary initial condition for the elliptic pde - if possible one that satisfies its boundary conditions.

If pdepe accepts it, it won't influence the overall solution in the sequel.

[Sharon]

Hi Torsten,

Thank you so much for answering my question! I give elliptic equation an initial condition which satisfies the boundary condition. But Matlab gives an error saying that "Failure at t=4.820300e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (8.881784e-16) at time t".

clc
clc
h=0.01;%space step size,
L=1;
T=10;
tau=0.001;%time step size
M=L/h;
x=linspace(0,L,M+1);
N=T/tau;
t=linspace(0,T,N+1);
m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);

Warning: Failure at t=6.785165e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t.

Warning: Time integration has failed. Solution is available at requested time points up to t=6.780000e-01.

u1 = sol(:,:,1);
u2 = sol(:,:,2);
plot(x,[u1(end,:);u2(end,:)]);

%%%
function [c,f,s] = pdefun(x,t,u,dudx)
chi=12;
c = [1; 0];
f = [(1*dudx(1))-chi*(u(1)/u(2))*(dudx(2)); dudx(2)];
s = [u(1)*(1-u(1)); -u(2)+u(1)];
end
%%
    function u0 = pdeic(x)
        u0 = [1+0.5*cos(pi*x); 1+0.5*cos(pi*x)];
    end
%%
    function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
        pl = [0; 0];
        ql = [1; 1];
        pr = [0; 0];
        qr = [1; 1];
    end

Could you please help me look at it? Thanks!

[Torsten]

Could you include your equations and boundary conditions in a mathematical notation so that it's possible to compare them with your settings in "pdepe" ?

[Sharon]

The equations are

The boundary conditions are

Thanks a lot for your help!

[Torsten]

I used the "correct" initial condition for v at t=0.

I thought "pdepe" should be able to adjust it automatically, but this is not the case.

syms x v(x)
u = 1+0.5*cos(pi*x);
eqn = diff(v,x,2)-v+u==0;
Dv = diff(v,x);
v = dsolve(eqn,[Dv(0)==0,Dv(1)==0])

v = 

v = matlabFunction(v);
u = matlabFunction(u);
%
h=0.01;%space step size,
L=1;
T=10;
tau=0.001;%time step size
M=L/h;
x=linspace(0,L,M+1);
N=T/tau;
t=linspace(0,T,N+1);
m = 0;
sol = pdepe(m,@pdefun,@(x)pdeic(x,u,v),@pdebc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
plot(x,[u1(end,:);u2(end,:)]);

%%%
function [c,f,s] = pdefun(x,t,u,dudx)
  chi=12;
  c = [1; 0];
  f = [dudx(1)-chi*u(1)/u(2)*dudx(2); dudx(2)];
  s = [u(1)*(1-u(1)); -u(2)+u(1)];
end
%%
function u0 = pdeic(x,u,v)
  u0 = [u(x); v(x)];
end
%%
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
  pl = [0; 0];
  ql = [1; 1];
  pr = [0; 0];
  qr = [1; 1];
end

[Sharon]

Thank you so much for modifying the code.

I am wondering how you get the "correct" initial value for v. You first solve with boundary condition to get v(x). Then you use v(x) as the initial value of and as the initial value of parabolic equation to solve the parabolic-elliptic system using pdepe. Is my understanding correct? Thanks!

[Torsten]

Yes.

The elliptic pde for v ( v''(x,0)-v(x,0)+u(x,0) = 0, v'(0,0)=v'(1,0) = 0 ) and the initial condition for u ( u(x,0) = 1+0.5*cos(pi*x) ) determine the initial condition for v at t=0 ( v(x,0) = 1 + 0.5*cos(pi*x)/(1+pi^2) ).

If your elliptic pde becomes more complicated such that an analytical solution for v(t=0) is no longer possible using the symbolic approach, you can also try to calculate v(t=0) using "bvp4c".

As said, "pdepe" usually manages to adjust an approximation to v(t=0) automatically, but seems to have difficulties in your case.

[Sharon]

Thank you a lot for such clear answers.

I am trying to use "bvp4c". But v at t=0 obtained by "bvp4c is some points, not a function. How to use such v as an initial function in pdepe? Thanks a lot for your help!

[Torsten]

You get a vector of values (x,v) from "bvp4c".

Define the function

vfun = @(xq) interp1(x,v,xq)

and pass it to "pdeic" as I did for the functions u and v in the symbolic approach.

Then - for a given x-value in "pdeic" - return

u0 = [1+0.5*cos(pi*x); vfun(x)];

to "pdepe".

[Torsten]

h=0.01;%space step size,
L=1;
T=10;
tau=0.001;%time step size
M=L/h;
x=linspace(0,L,M+1);
N=T/tau;
t=linspace(0,T,N+1);
%
u = @(x)1+0.5*cos(pi*x);
ux = @(x)-0.5*pi*sin(pi*x);
bvpfcn = @(x,v)[v(2);v(1)-u(x)];
bcfcn = @(va,vb)[va(2);vb(2)];
guess = @(x)[u(x);ux(x)];
solinit = bvpinit(x,guess);
sol = bvp4c(bvpfcn, bcfcn, solinit);
%plot(sol.x,sol.y(1,:))
v = @(xq)interp1(sol.x,sol.y(1,:),xq);
%
m = 0;
sol = pdepe(m,@pdefun,@(x)pdeic(x,u,v),@pdebc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
plot(x,[u1(end,:);u2(end,:)]);

%%%
function [c,f,s] = pdefun(x,t,u,dudx)
  chi=12;
  c = [1; 0];
  f = [dudx(1)-chi*u(1)/u(2)*dudx(2); dudx(2)];
  s = [u(1)*(1-u(1)); -u(2)+u(1)];
end
%%
function u0 = pdeic(x,u,v)
  u0 = [u(x); v(x)];
end
%%
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
  pl = [0; 0];
  ql = [1; 1];
  pr = [0; 0];
  qr = [1; 1];
end

[Sharon]

Thank you so much for your help!

But the graph got doesn't match the theoretical result. Theoretically, the solution for u should near 1 and behaves like cosin function. I am thinking maybe my equations cannot be solved in this way. Is it possible to solve it using PDE toolbox?

[Torsten]

As far as I know, the spatial dimension for the PDE toolbox is at least 2. So I think you cannot use it for your problem.

This is in parts confirmed by the following contribution:

https://de.mathworks.com/matlabcentral/answers/91863-can-i-solve-1-dimensional-problems-using-the-partial-differential-equation-pde-toolbox

Do you have a link to the "theoretical result that the solution should be near 1 and behave like a cosinus function" ?

[Sharon]

Yes. I see it. PDE toolbox may not work for my case.

I don't have a link. Actually, this result has not been proved rigorously. We are working on it. This is our conjecture. We believe this conjecture is true. I think maybe this is pdede solver's problem that cannot deal with such case?

[Torsten]

I think maybe this is pdede solver's problem that cannot deal with such case.

Maybe. But I doubt that a different solver will perform better.

Your system seems to be a challenge because reformulating the problem for pdepe gives a very different result.

Differentiating d/dx (du1/dx - chi*u1/u2*du2/dx) and inserting d^2u2/dx^2 = u2 - u1 leads to the following code with the following solution. Maybe there is an error in my reformulation, but I can't find any.

syms x v(x)
u = 1+0.5*cos(pi*x);
eqn = diff(v,x,2)-v+u==0;
Dv = diff(v,x);
v = dsolve(eqn,[Dv(0)==0,Dv(1)==0]);
v = matlabFunction(v);
u = matlabFunction(u);
%
h=0.01;%space step size,
L=1;
T=10;
tau=0.001;%time step size
M=L/h;
x=linspace(0,L,M+1);
N=T/tau;
t=linspace(0,T,N+1);
m = 0;
sol = pdepe(m,@pdefun,@(x)pdeic(x,u,v),@pdebc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
plot(x,[u1(end,:);u2(end,:)])

function [c,f,s] = pdefun(x,t,u,dudx)
chi = 12;
c = [1; 0];
f = [dudx(1);dudx(2)];
s = [-chi*((dudx(1)*u(2)-u(1)*dudx(2))/u(2)^2 * dudx(2) + u(1)/u(2)*(u(2)-u(1)))+u(1)*(1-u(1));-u(2)+u(1)];
%f = [dudx(1)-chi*u(1)/u(2)*dudx(2); dudx(2)];
%s = [u(1)*(1-u(1)); -u(2)+u(1)];
end
%%
function u0 = pdeic(x,u,v)
u0 = [u(x); v(x)];
end
%%
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
pl = [0; 0];
ql = [1; 1];
pr = [0; 0];
qr = [1; 1];
end

[Sharon]

Thank you so much for your effort. I checked it and also didn't find any error in your reformulation. Is it because the "correct" initial function is not correct?

[Torsten]

Is it because the "correct" initial function is not correct?

No. The initial condition for v is correct.

pdepe just gets different results for the two equivalent formulations of the problem. But why ?

I think the reason is the boundary condition.

It's not that simple in the first formulation for the solver to detect that setting

f = [dudx(1)-chi*u(1)/u(2)*dudx(2); dudx(2)]
ql(1) = qr(1) = ql(2) = qr(2) = 1

really leads to

du/dx = dv/dx = 0

on both ends.

You defined

 f = [dudx(1)-chi*u(1)/u(2)*dudx(2); dudx(2)]

Together with

ql(1) = qr(1) = ql(2) = qr(2) = 1

this gives

du/dx - chi*u/v*dv/dx = 0 and dv/dx = 0

at both ends since the boundary condition is

p+q*f = 0.

This is a linear system in du/dx and dv/dx which is solved by du/dx = dv/dx = 0.

But I don't know how the pdepe solver internally works. Judging from the different result, I guess that the two boundary condition settings are interpreted differently.

[Sharon]

It makes sense. Thank you so much!