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
2 个评论
Torsten
2019-1-16
Start with first-order upwind.
Malak Galal
2019-1-18
编辑：Malak Galal
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
2019-1-24
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
Torsten
2019-1-31
> 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
2019-1-31
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.
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
2019-1-31
Torsten
2019-1-31
Correct.
Malak Galal
2019-1-31
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
更多回答（3 个）
Malak Galal
2019-2-22
0 个投票
1 个评论
Torsten
2019-2-22
Hello Malak,
nx and nt can be chosen completely independent from each other.
Best wishes
Torsten.
类别
在 帮助中心 和 File Exchange 中查找有关 Linear Algebra 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
