%clear all; %close all; % Program mtc (method of characteristic) % % This program solves the initial value problem % % u_t + uu_x = 0, u(x,0) = f(x) % % using the method of characteristics. There are four choices of % initial data. Choices 1 and 2 should develop a shock, but instead % 'roll over' into a curve which no longer defines a function. % Choice 3 is a rarefaction wave. The initial data for choice 4 % is a small perturbation of the constant function f(x) = 1. In this % case the program computes the solution of the nonlinear equation % and the solution of the linearized problem which is u_t + u_x = 0. % At run time user enters the choice of data, 1,2,3, or 4, and % the times t1,t2,t3,and t4 at which the solution is to be viewed. % The solution profiles at times t = 0,t1,t2,t3,t4 are all plotted on % the interval [-1 6]. The solution profiles of both linear and nonlinear % problems are plotted together in the case of data option 4. disp(' ') disp('Enter 1 for decreasing profile ') disp('Enter 2 for another decreasing profile ') disp('Enter 3 for increasing profile ') disp('Enter 4 for the small perturbation of f(x)= 1 ') disp('Enter 5 for po') m = input( 'enter the choice of initial data: 1,2,3,4 or 5 ') t = input(' Enter the times in the form [t1 t2 t3 t4 ] ') x = -5:.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); elseif m == 2 ymiddle = .25*(x+1).*(x-2).^2; y = (x < 0) + ( (x < 2.001 ) - (x < 0) ).*ymiddle ; elseif m == 3 ymiddle = .5 - (x+1).^2 .*(x-2)/8 ; y = .5*(x < -1) + ( ( x < 1.001 ) - ( x < -1) ).*ymiddle + (x > 1); elseif m == 4 ymiddle = .05*x.^3 - .15*x.^2 + 1.1; y = 1.1*(x < 0) + ( (x < 2.001) - (x < 0) ) .*ymiddle + .9*(x > 2); elseif m==5 %y=pi/2+atan(3*x); %y=pi/2-atan(x); %ymiddle=1-0.5*x; %y = (x < 0)*1 + ( (x < 1.001 ) - (x< 0) ).*ymiddle + .5*(x > 1); ymiddle=1-.5*x; y = (x < 0)*1 + ( (x < 1.001 ) - (x< 0) ).*ymiddle + 0.5*(x > 1); end high = max(y); low = min(y); gap =.25*(high - low); x1 = x + t(1)*y; x2 = x + t(2)*y; x3 = x + t(3)*y; x4 = x + t(4)*y; figure plot(x,y, 'r') axis([-5,6, low-gap, high+gap]) hold on plot(x1,y,'b-o') plot(x2,y,'y-+') plot(x3,y,'k-.') plot(x4,y,'g--') hold off say1=['t=' num2str(t(1))]; say2=['t=' num2str(t(2))]; say3=['t=' num2str(t(3))]; say4=['t=' num2str(t(4))]; legend('t=0',say1,say2,say3,say4,'Location', 'Southeast' ) if m == 4 xl1 = x + t(1); xl2 = x + t(2); xl3 = x + t(3); xl4 = x + t(4); plot(xl1,y,'g') plot(xl2,y,'g') plot(xl3,y,'g') plot(xl4,y,'g') else end %set(gca, 'XLim',[-1 10]) hold off %legend('1','2','3','4','5')