solving nonlinear wave equation
6 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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
0 Commenti
Risposta accettata
Torsten
il 12 Giu 2024
Modificato: Torsten
il 13 Giu 2024
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 Commenti
Torsten
il 13 Giu 2024
Modificato: Torsten
il 13 Giu 2024
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)
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!