First order partial differential equations system - Numerical solution

60 ビュー (過去 30 日間)
Malak Galal
Malak Galal 2019 年 1 月 16 日
コメント済み: Sreejath S 2020 年 10 月 26 日
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
Torsten 2019 年 1 月 16 日
Start with first-order upwind.
Malak Galal
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

サインインしてコメントする。

採用された回答

Torsten
Torsten 2019 年 1 月 18 日
編集済み: 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 件のコメント
Malak Galal
Malak Galal 2019 年 1 月 31 日
Okay. Thank you very much for your help!
Sreejath S
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

サインインしてコメントする。

その他の回答 (2 件)

Malak Galal
Malak Galal 2019 年 1 月 21 日
Hi Torsten,
Thank you very much for your help!
Malak

Malak Galal
Malak Galal 2019 年 2 月 22 日
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 件のコメント
Torsten
Torsten 2019 年 2 月 22 日
Hello Malak,
nx and nt can be chosen completely independent from each other.
Best wishes
Torsten.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by