%close all; %clear all; % Matlab program cl (conservation law) % % This program uses an upwind finite difference code to solve the % initial value problem u_t + (F(u))_x= 0, u(x,0) = f(x) . % The u_t is approximated by a forward difference and (F(u))_x % by a backward difference. The method is dissipative and % tend to smear out the shocks. % User must provide a flux function F(u). It is required that % c(u) = F'(u) > 0 for all values of u in the problem. As an % example F(u) = .5*u^2 is given. % It is further assumed that the left boundary value is constant. % The program calls the function "stepcl" which integrates the % finite difference scheme "nsteps" time steps foward in time. % There are six choices of initial data. The first two are the same % as for the program mtc. They % give rise to shock waves. The third and fourth are simple step % functions. The fifth and sixth are % combinations of humps. % User may make his own choice of initial data in the % function f.m. The function must be array smart. % At run time user must make the choice of data: 1,2,3,4,5,6, or 7. % User must also enter the size of the spatial step size delx, and the % size of the time step delt. delt and delx must be chosen to % satisfy the CFL condition. % Also at run time the user must enter the number n1 of time steps to % time t1, the number of time steps n2 from time t1 to time t2, etc. % For example if delt = .05, then n1 = n2 = 20, and n3 = 10, n4 = 20 % yields the times t1 = 1, t2 = 2, t3 = 2.5, and t4 = 3.5. % % The program saves for viewing snapshots of the solution at times % t = 0 (snap0), t = t1 (snap1), t = t2 (snap2), etc. All five snapshots % are plotted together on the interval -1 < x < 6. To plot one of % them separately, say snap2, use the command % % plot(x,snap2) % % To see two of them together, say snap2 and snap4, use the command % % plot(x,snap2,x,snap4) % disp(' ') disp('Enter 1 for decreasing profile - develops shock ') disp('Enter 2 for another of the same ' ) disp('Enter 3 for step shock wave ') disp('Enter 4 for centered rarefaction wave ') disp('Enter 5 for single hump ') disp('Enter 6 for hump and valley ') disp('Enter 7 for your choice ') m = input('Enter 1,2,3,4,5,6, or 7 ') delx = input('Enter the spatial step delta x ') delt = input('Enter the time step delta t ') N = input(' Enter the number of time steps in the form [n1, n2, n3, n4] ') n1 = N(1); t1 = n1*delt n2 = N(2); t2 = (n1+n2)*delt n3 = N(3); t3 = (n1 + n2 +n3)*delt n4 = N(4); t4 = (n1 + n2 + n3 + n4)*delt rho = delx/delt; x = -1:delx:6; jmax = length(x); if m == 1 umiddle = 1.0 - .125*x.^2.*(3.0 - x); u = (x < 0) + ( (x < 2.001 ) - (x< 0) ).*umiddle + .5*(x > 2); elseif m == 2 umiddle = .25*(x+1).*(x-2).^2; u = (x < 0) + ( (x < 2.001 ) - (x < 0) ).*umiddle ; elseif m == 3 u = .99999*(x< delx/2); elseif m == 4 u = .99999*(x > -delx/2); elseif m == 5 u = .5*exp(-2.0*(x-1).^2) + .25; elseif m == 6 u = 2.0*(x-1).*exp(-4.0*(x-1).^2) + .5; elseif m == 7 u = f(x); end snap0 = u; j = [2:jmax]; for n =1:n1 u(j) = u(j) - (flux(u(j)) - flux(u(j-1)) )/rho; end snap1 = u; for n = 1:n2 u(j) = u(j) -( flux(u(j)) - flux(u(j-1)) )/rho; end snap2 = u; for n = 1:n3 u(j) = u(j) -( flux(u(j)) - flux(u(j-1)) )/rho; end snap3 = u; for n = 1:n4 u(j) = u(j) -( flux(u(j)) - flux(u(j-1)) )/rho; end snap4 = u; %figure plot(x, snap0,'k') %break; hold on plot(x,snap1,'c') hold on plot(x,snap2,'y') hold on plot(x,snap3,'b') hold on plot(x,snap4,'m') %hold off legend('t0', 't1', 't2', 't3', 't4', 't5','NorthWestOutside')