First order partial differential equations system - Numerical solution

Hello everyone!

I would like to solve a first order partial differential equations (2 coupled equations) system numerically. I just want to make sure that my thoughts are correct.

So firstly, I will start by doing a discretization to each of the two equations and then I will use ode15s to solve the ordinary differential equations that I got from the first step. Am I right? Which method of discretization do you recommend me to use?

Note: My equations are similar to the linear advection equation but with a source term.

Thank you very much.

Malak

[Torsten]

Start with first-order upwind.

[Malak Galal]

Thank you very much for your reply.

Is there any document or website that could be of use to me?

I am a bit lost because I have been trying to understand how to use the upwind method to solve my system, but I couldn't really find anything similar to my system. I saw some examples about the advection equation, but they all deal with one equation and with no source term.

My problem is that the source terms in my system couple the two equations.

I would appreciate your help very much! Thanks in advance.

Malak

Answer 1

[Torsten]

%program solves
% dy1/dt + v1*dy1/dx = s1*y1*y2
% dy2/dt + v2*dy2/dx = s2*y1*y2
% y1(0,x) = y2(0,x) = 0
% y1(t,1) = y2(t,0) = 1
% for 0 <= t <= 400

function main
  nx = 500;
  y1 = zeros(nx,1);
  y1(end) = 1.0;
  y2 = zeros(nx,1);
  y2(1) = 1.0;
  ystart = [y1;y2];
  tstart = 0.0;
  tend = 400;
  nt = 41;
  tspan = linspace(tstart,tend,nt);
  xstart = 0.0;
  xend = 1.0;
  x = linspace(xstart,xend,nx);
  x = x.';
  v1 = -0.005;
  v2 = 0.0025;
  v = [v1 v2];
  s1 = 3.0e-3;
  s2 = -4.0e-3;
  s = [s1 s2];
  [T,Y] = ode15s(@(t,y)fun(t,y,nx,x,v,s),tspan,ystart);
  plot(x,Y(20,1:nx), x,Y(20,nx+1:2*nx))
end

function dy = fun(t,y,nx,x,v,s)
  y1 = y(1:nx);
  y2 = y(nx+1:2*nx);
  dy1 = zeros(nx,1);
  dy2 = zeros(nx,1);
  dy1(1:nx-1) = -v(1)*(y1(2:nx)-y1(1:nx-1))./(x(2:nx)-x(1:nx-1)) + s(1)*y1(1:nx-1).*y2(1:nx-1);
  dy1(nx) = 0.0;
  dy2(1) = 0.0;
  dy2(2:nx) = -v(2)*(y2(2:nx)-y2(1:nx-1))./(x(2:nx)-x(1:nx-1)) + s(2)*y1(2:nx).*y2(2:nx) ;
  dy = [dy1;dy2];
end

[Malak Galal]

Hi Torsten,

I would like to thank you again very much for your help. I was working on understanding exactly how the code works to be able to solve my equations; so of course some questions popped up.

1) The step where you put both column vectors (y1 and y2) in one column vector (y); is it because the ode function cannot take more than one column vector?

2) When I reverse the boundary conditions, for example like having y1(1)=1 and y2(end)=0, why doesn't the code give proper solutions anymore? It seems like the code is written in such a way that y1 has to propagate from right to left and y2 from left to right, and if we reverse them, it won't work. Could you please clarify this matter to me?

Thanks a lot!

Malak

[Torsten]

ad 1)

Yes, the MATLAB ODE integrators require one single column vector for the solution vector (y) and its time derivatives (dy). Thus if you have to solve several PDEs, you will have to arrange the variables somehow in one vector.

ad 2)

Advection equations only permit one boundary condition, not two.

If v < 0, the solution propagates from right to left. In this case, a boundary condition has to be set at the right end of the solution domain and the derivative dy/dx in point x(i) should be approximated as

(dy/dx)(xi) ~ (y(x(i+1))-y(x(i))/(x(i+1)-x(i)) 

to get a stable 1st-order accurate solution. As you can see, no boundary condition in the left endpoint is necessary because an ODE can be written for dy(1)/dt.

If v > 0, the solution propagates from left to right. In this case, a boundary condition has to be set at the left end of the solution domain and the derivative dy/dx in point x(i) should be approximated as

(dy/dx)(xi) ~ (y(x(i))-y(x(i-1))/(x(i)-x(i-1)) 

to get a stable 1st-order accurate solution. As you can see, no boundary condition in the right endpoint is necessary because an ODE can be written for dy(end)/dt.

Best wishes

Torsten.

[Malak Galal]

I see. Thanks a lot for your reply.

If I want to test a different input, for instance instead of having y2(1)= 1 we would have y2(2:100)= 1, which is basically a pulse; why does ODE15s take so much time to run? I tried using ODE45 instead, but I know it should be used with non-stiff equations. It, however, runs much faster, but keeps smoothening the pulse and for some reason the pulse keeps decaying even without inducing any coupling.

Do you know why would the pulse decay even without any coupling? Why would the pulse change its shape with time? Is it because of the use of ODE45?

Thanks in advance!

Malak

[Torsten]

> If I want to test a different input, for instance instead of having y2(1)= 1 we would have y2(2:100)= 1, which is basically a pulse;

I don't understand what you try to set here. Do you want to change the initial profile of y2 such that it is 1 for 1 <=i <= 100 and 0 for 101 <=i <=500 ?

Best wishes

Torsten.

[Malak Galal]

Yes, exactly. I am sorry if I was unclear.

[Torsten]

And why do you talk about a "pulse" for this ? Usually a pulse is something that appears for a short instant and disappears immediately afterwards.

Furthermore: Which boundary condition do you want to set at x=xstart for y2 ?

[Malak Galal]

I am talking about a pulse that has a specific width, not just a dirac. What I want to achieve is a propagating pulse. This pulse should propagate from 0 to end in space and accordingly should have a different position at each instance of time.

My boundary condition would be xstart=0 as it already is. So instead of y2(1)=1 which is only one point in space equal to 1, I want to have more than one point equal to 1 and the rest is 0.

I already did that by doing y2(2:100)=1; it gives the shape of a rectangular pulse with a specific width. At different time instances, this pulse starts to decay and changes its shape to look like a Gaussian pulse as shown on the figure below. The y-axis is the amplitude that I set to 1, and the x-axis are the 500 points in space. And this is how the pulse evolutes with time.

I have absolutely no idea why it decays and transforms into a Gaussian pulse, even though I disabled the coupling. It should just be the same pulse at all time instances.

[Torsten]

function main
  nx = 20000;
  y = zeros(nx,1);
  y(2:4000) = 1.0;
  ystart = y;
  tstart = 0.0;
  tend = 400;
  nt = 41;
  tspan = linspace(tstart,tend,nt);
  xstart = 0.0;
  xend = 1.0;
  x = linspace(xstart,xend,nx);
  x = x.';
  v = 0.0025;
  [T,Y] = ode45(@(t,y)fun(t,y,nx,x,v),tspan,ystart);
  plot(x,Y(10,:),x,Y(20,:),x,Y(30,:))
end



function dy = fun(t,y,nx,x,v)
  dy = zeros(nx,1);
  dy(1) = 0.0;
  dy(2:nx) = -v*(y(2:nx)-y(1:nx-1))./(x(2:nx)-x(1:nx-1)) ;
end

Better ?

That's why I said: Start with first-order upwind.

If you want sharp profiles with less grid points, you will have to discretize the first derivative term with a second-order upwind method (Lax-Wendroff etc.).

[Malak Galal]

Okay! So, basically it has nothing to do with ODE45, it all depends on the discretization method; am I correct?

[Torsten]

Correct.

[Malak Galal]

Okay. Thank you very much for your help!

[Sreejath S]

Hi Torsten,

Sorry for commenting on this post, rather than asking my own question. I have a very similar system of coupled first order hyperbolic PDES. But my system has a major diference from the example code you posted.

The first equation has dy1/dt and dy2/dx, and the second equation has dy2/dt and dy1/ds. Kindly note the change (marked bold) in the following lines from the original problem you posted:

% dy1/dt + v1*dy2/dx = s1*y1*y2

% dy2/dt + v2*dy1/dx = s2*y1*y2

% y1(0,x) = y2(0,x) = 0

% y1(t,1) = y2(t,0) = 1

% for 0 <= t <= 400

Can this type of coupled problem be solved using MOL?

I ran the code with this change, but I could not get a solution. The solution is not converging.Kindly help in this regard.

Sincerly,

Sreejath

Answer 2

[Malak Galal]

Hi Torsten,

Thank you very much for your help!

Malak

Answer 3

[Malak Galal]

Hi Torsten,

I have a question concerning the number of points in space and time (nx and nt). Using a proper discretization method, if I want to use for instance nx=1e5, does nt have to be the same as nx or close to it? It seems to me that when I increase the number of points in space, the number of points in time have to be increased as well. Could you please help me understand the relationship between nx and nt?

Thank you very much!

Malak

[Torsten]

Hello Malak,

nx and nt can be chosen completely independent from each other.

Best wishes

Torsten.