solving nonlinear wave equation
显示 更早的评论
I need to solve this PDE below.
This is quite similar to the Burgers' equation. However, I can't understand how to get the exact solution to this equation. Please give an explicit solution to a given initial condition. Any initial condition is fine, just as long as the explicit solution is given. I will use it to check if my FDM code is correct.
In addition, if anyone has the code to solving this equation numerically, please post it here.
Thanks in advance!
below is my code.
for i = 2:(step(1)+1) % t
for j = 2:(step(2)+1) % x
rho(i, j) = rho(i-1, j) - v_max*(1-2*(rho(i-1, j)/rho_max))*(dt/dx)*(rho(i-1, j) - rho(i-1, j-1));
end
end
%----------- set up ---------------
tspan = [0, 3]; % t
xspan = [-10, 10]; % x
step = [1000, 1000]; % step(1) : t, step(2) : x
t_values = linspace(tspan(1), tspan(2), step(1)+1);
x_values = linspace(xspan(1), xspan(2), step(2)+1);
v_max = 5; rho_max = 5; % vmax, rhomax
dt = (tspan(2) - tspan(1)) / step(1);
dx = (xspan(2) - xspan(1)) / step(2);
rho = zeros(step(1)+1, step(2)+1);
rhoo = zeros(step(1)+1, step(2)+1);
for j = 1:(step(1)+1)
rho(1, j) = 1/(1 + exp(j*dt));
end
%------------- (FDM) -------------
for i = 2:(step(1)+1) % t
for j = 2:(step(2)+1) % x
rho(i, j) = rho(i-1, j) - v_max*(1-2*(rho(i-1, j)/rho_max))*(dt/dx)*(rho(i-1, j) - rho(i-1, j-1));
end
end
% for i = 1:(step(1)+1)
% for j = 1:(step(2)+1)
% rhoo(i, j) =
% end
% end
%------------ animation --------------
filename = 'animation.gif';
figure;
for i = 1:(step(1)+1)
plot(x_values, rho(i, :));
xlabel('Position');
ylabel('Density');
title(['Time = ', num2str(t_values(i))]);
drawnow;
frame = getframe(gcf);
img = frame2im(frame);
[imind, cm] = rgb2ind(img, 256);
if i == 1
imwrite(imind, cm, filename, 'gif', 'Loopcount', inf, 'DelayTime', dt);
else
imwrite(imind, cm, filename, 'gif', 'WriteMode', 'append', 'DelayTime', dt);
end
end
采纳的回答
Using the method of characteristics, you get the equations
dt/ds = 1
dx/ds = v_max*(1-2*rho/rho_max)
drho/ds = 0
You should be able to solve these 3 ordinary differential equations analytically and get the analytical solution for your PDE. This is easy for your initial conditions for rho since rho(t=0,x) is monotonically decreasing in x. That guarantees that the characteristic lines cannot cross.
A first indication for the precision of your code is the line that separates the regions rho = 0 and rho > 0. With your settings, it should be t = 0.25*(x+10). Here is the code to check this. I also included a plot of the maximum rho over time which should remain rho = 0.5 throughout:
%----------- set up ---------------
tspan = [0, 3]; % t
xspan = [-10, 10]; % x
step = [1000, 1000]; % step(1) : t, step(2) : x
t_values = linspace(tspan(1), tspan(2), step(1)+1);
x_values = linspace(xspan(1), xspan(2), step(2)+1);
v_max = 5; rho_max = 5; % vmax, rhomax
dt = (tspan(2) - tspan(1)) / step(1);
dx = (xspan(2) - xspan(1)) / step(2);
rho = zeros(step(1)+1, step(2)+1);
for j = 1:(step(1)+1)
rho(1, j) = 1/(1 + exp(0.15*(x_values(j)+10)));
end
for i = 2:(step(1)+1) % t
for j = 2:(step(2)+1) % x
rho(i, j) = rho(i-1, j) - v_max*(1-2*(rho(i-1, j)/rho_max))*(dt/dx)*(rho(i-1, j) - rho(i-1, j-1));
end
end
for i = 1:step(1)+1
[rhomax(i),idx] = max(rho(i,:));
x(i) = x_values(idx);
end
figure(1)
hold on
plot(t_values,4*t_values-10,'r--')
plot(t_values,x,'b') % Should be t = 0.25*(x+10) (moving front)
hold off
grid on
x(end)
figure(2)
hold on
plot(t_values,0.5*ones(size(t_values)),'r--')
plot(t_values,rhomax,'b')
hold off
grid on
For t = 3, x should be equal to 3/0.25 - 10 = 2. Here, it is approximately x = 3.62 .
2 个评论
Student
2024-6-13
The characteristics starting in (-10,t) for t>0 and in (x,0) for x > 0 intersect. So information about the value for rho in (x,t) comes from two sides. I'm quite sure that the line where the maximum of rho occurs is t = 0.25*(x+10), but the value will most probably decrease from 0.5 because of dissipation.
If you want a reliable result, you should apply CLAWPACK to the rewritten equation
drho/dt + df(rho)/dx = 0
with
f(rho) = v_max*rho*(1 - rho/rho_max)
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 PDE Solvers 的更多信息
Translated by 
