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

Answer 1

[Torsten]

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)

ans = 3.6200

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 .

[Student]

Thanks, it helped me a lot!!

[Torsten]

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)

https://www.clawpack.org/