% Matlab program shocks % % This program solves the initial value problem % % u_t + uu_x = 0, u(x,0) = f(x) % % using the method of characteristics, modified to capture shocks when % they form. There are four choices of initial data. Choices 1 and 2 % develop a shock. Choice 3 is a step function which also yields a % shock. % At run time user enters the choice of data, 1,2, or 4 and % the times t1,t2,t3,and t4 at which the solution is to be viewed. % Solutions are all plotted on the interval -1 < x < 6. disp( ' ') disp(' Enter 1 for decreasing profile - develops shock ') disp(' Enter 2 for another decreasing profile ') disp(' Enter 3 for step ') m = input( 'enter the choice of initial data: 1,2, or 3 ') t = input(' Enter the times in the form [t1 t2 t3 t4] ') t1 = t(1); t2 = t(2); t3 = t(3); t4 = t(4); x = -1:.05:6; n = length(x); if m == 1 % ymiddle = 1.0 - .125*x.^2.*(3.0 - x); %y = (x < 0) + ( (x < 2.001 ) - (x< 0) ).*ymiddle + .5*(x > 2); ymiddle=1-.5*x; y = (x < 0)*1 + ( (x < 1.001 ) - (x< 0) ).*ymiddle + 0.5*(x > 1); elseif m == 2 ymiddle = .25*(x+1).*(x-2).^2; y = (x < 0) + ( (x < 2.001 ) - (x < 0) ).*ymiddle ; elseif m == 3 y = (x < 0.0001)+(x>0)*0.5; end if m == 3 x1 = x + 0.75*t1; x2 = x + 0.75*t2; x3 = x + 0.75*t3; x4 = x + 0.75*t4; else s = y(41); x1 = zeros(size(x)); for j = 1:141 r = (x(j) + t1*y(j)-(1 +s*t1))*(x(j)-1); if r > 0 x1(j) = x(j) + t1*y(j); else x1(j) = 1+s*t1; end end x2 = zeros(size(x)); for j = 1:141 r = (x(j)+ t2*y(j) - (1+s*t2))*(x(j)-1); if r > 0 x2(j) = x(j) +t2*y(j); else x2(j) = 1 + s*t2; end end x3 = zeros(size(x)); for j = 1:141 r = (x(j) + t3*y(j)-(1+s*t3))*(x(j)-1); if r > 0 x3(j) = x(j) +t3*y(j); else x3(j) = 1 + s*t3; end end x4 = zeros(size(x)); for j = 1:141 r = (x(j) + t4*y(j)-(1+s*t4))*(x(j) -1); if r > 0 x4(j) = x(j) + t4*y(j); else x4(j) = 1 +s*t4; end end end high = max(y); low = min(y); gap = .25*(high - low); plot(x,y,'r') axis([-1, 6, low-gap, high+gap]) hold on plot(x1,y,'g') plot([-1, x1(21)], [y(1), y(1)], 'g') plot(x2,y,'g') plot([-1 x2(21)], [y(1), y(1)], 'g') plot(x3,y,'g') plot([-1, x3(21)], [y(1),y(1)],'g') plot(x4,y,'g') plot([-1 x4(21)], [y(1),y(1)],'g') hold off