First order partial differential equations system - Numerical solution

Question

Malak Galal 2019-1-16

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/440015-first-order-partial-differential-equations-system-numerical-solution

编辑： Torsten 2025-1-18

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

2 个评论
显示无隐藏无

Torsten 2019-1-16

Start with first-order upwind.

Malak Galal 2019-1-18

编辑：Malak Galal 2019-1-18

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 2019-1-18

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/440015-first-order-partial-differential-equations-system-numerical-solution#answer_356975

编辑：Torsten 2019-1-18

%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

12 个评论
显示 10更早的评论隐藏 10更早的评论

Torsten 2019-1-25

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 2019-1-31

编辑：Malak Galal 2019-1-31

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.

Malak Galal 2019-1-31

Okay. Thank you very much for your help!

Sreejath S 2020-10-26

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 2019-1-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/440015-first-order-partial-differential-equations-system-numerical-solution#answer_357275

Hi Torsten,

Thank you very much for your help!

Malak

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 3

Malak Galal 2019-2-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/440015-first-order-partial-differential-equations-system-numerical-solution#answer_362289

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2019-2-22

Hello Malak,

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

Best wishes

Torsten.

请先登录，再进行评论。

Answer 4

Raghunath 2025-1-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/440015-first-order-partial-differential-equations-system-numerical-solution#answer_1557658

dy/dx+3x=0,y{0)=1

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2025-1-18

编辑：Torsten 2025-1-18

Do you know the antiderivatives y(x) of -3*x ? Can you adjust it such that y(0) = 1 ?

请先登录，再进行评论。

First order partial differential equations system - Numerical solution

2 个评论
显示无隐藏无

采纳的回答

12 个评论
显示 10更早的评论隐藏 10更早的评论

更多回答（3 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

First order partial differential equations system - Numerical solution

2 个评论 显示 无隐藏 无

采纳的回答

12 个评论 显示 10更早的评论隐藏 10更早的评论

更多回答（3 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

12 个评论
显示 10更早的评论隐藏 10更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论