% viscocmd.m % % Created by Victor Camacho on 5/20/06 % % In this script we consider the PDE in the following form (s = sigma) % % u_t = s_x - (u^2/2)_x % % s_t = (a^2+s)*u_x - b*s - u*s_x % % Where a and b are assumed to be constant. % % In this form it is easy to apply a finite difference approximation at % each time step for both dependent variables. % % This script serves as a series of callback routines from the parent file % visco.m which sets up the Graphical User Interface for the main program. % As an input this function takes the command, given by visco.m, and % excecutes the proper set of functions based on the command's value. function viscocmd(cmd) % These global structures are more rigorously defined in visco.m, however % several of them will be constructed in this program and definitions of % these structures' components will be described. global myfig; global sol; global sol2; global grid; global playopts; global parameters; global cancel; switch (cmd) % If the user has closed the program's figure, this command will be % sent and so we must delete the figure objects we've created and clear % all global variables from memory. case 'close' delete(myfig.fig1); clear global myfig sol grid playopts parameters; fig2 = findobj('Tag','analyzewindow'); if (~isempty(fig2)) delete(fig2); end fig3 = findobj('Tag','analyzewindow2'); if (~isempty(fig3)) delete(fig3); end return; % If the user only wishes to close the analyzing window, the program % continues to run and we must preserve the contents of the window. So % the window is simply set to invisible upon closure. case 'close2' fig2 = findobj('Tag','analyzewindow'); if (~isempty(fig2)) set(myfig.fig2,'visible','off'); end return; case 'close3' fig3 = findobj('Tag','analyzewindow2'); if (~isempty(fig3)) set(myfig.fig3,'visible','off'); end return; % After the PDE has been solved, the figure's "Play" button is enabled % and when clicked, the 'play' command is sent to this program. We % obtain the user's time window for playback and see that it makes % sense. If it does we run the play() function, described in detail % later on, with the appropriate time settings. Otherwise an error is % displayed in the status window. case 'play' cancel = false; start = str2num(get(myfig.text.playstart,'string')); finish = str2num(get(myfig.text.playend,'string')); if (isempty(start) | isempty(finish)) errmsg('Bad play time settings!!'); return; end play(start,finish); % After the user has entered all of the desired specifications for % simulating the PDE, they hit the button "Run" to run the simulation. % This sends the 'run' command to this file. We begin by checking for % errors in the data which the user has entered. If an error is % detected it is displayed in the status bar and the function returns. % Otherwise, we take all of the user input and store it in the % necessary global structures to be used in subsequent function. Next % the contents of the sol structure are deleted to make room for the % new solution. Next the initial conditions are calculated and placed % in the first columb of the solution matrices. If something went % wrong in these calculations an error is returned to the status bar. % Otherwise, the program calls simulateSolution() which calculates the % solution of the PDE given the initial conditions and stores the vital % information to the sol structure. Finally, the solution is % automatically played back once upon the completion of this % calculation and the "Play" button and "Continue" button are enabled. case 'run' cancel = false; if (~validconditions(myfig)) return; end settings(); clearSolution(); if (cancel == true) return; end ic = initialConditions(); if (isempty(ic.u)) errmsg('Invalid initial conditions!!'); return; end method1 = get(myfig.menu,'Value'); method2 = get(myfig.menu2,'Value')-1; sol = simulateSolution(ic,method1); if (method2 ~= 0 & method2 ~= method1) sol2 = simulateSolution(ic,method2); end if (cancel == true) return; end play(0,grid.tfinal); activateButtons(); % After a simulation has been run, the user can continue solving the % PDE for subsequent times using the state at the final time as the % intitial condition. When the user hits the "Continue" button then % the 'continue' command is passed to this program. We first check % that the user entered a valid continue time. If so then we set the % initial conditions appropriately and simulate the solution as before % only with a new set of initial conditions. Then this new solution is % appended to the old one and stored in the sol structure. The mesh on % which the solution is defined changes as well so it is modified % appropriately by changing the grid structure. Finally, the solution % is played back starting from the previous ending time and ending at % the updated ending time. case 'continue' cancel = false; try time = eval(get(myfig.text.continue,'string')); catch errmsg('Invalid continue time'); return; end grid.T = floor(time*grid.tprec); ic = initialConditions(); method1 = get(myfig.menu,'Value'); method2 = get(myfig.menu2,'Value')-1; s = simulateSolution(ic,method1); if (method2 ~= 0 & method2 ~= method1) s2 = simulateSolution(ic,method2); sol2.u = [sol2.u s2.u]; sol2.s = [sol2.s s2.s]; if (parameters.wave) sol2.uw = [sol2.uw s2.uw]; sol2.sw = [sol2.sw s2.sw]; end end sol.u = [sol.u s.u]; sol.s = [sol.s s.s]; if (parameters.wave) sol.uw = [sol.uw s.uw]; sol.sw = [sol.sw s.sw]; end [trash grid.T] = size(sol.u); clear s; clear trash; grid.tfinal = grid.tfinal+time; set(myfig.text.tfinal,'string',num2str(grid.tfinal)); play(grid.tfinal-time,grid.tfinal); set(myfig.text.playstart,'string','0'); set(myfig.text.playend,'string',num2str(grid.tfinal)); return; % One of the checkboxes in the figure specifies whether the user would % prefer to define the initial conditions as a wave perturbation or % not. When this checkbox is clicked (turned on or off) this function % is sent the command 'pert'. All we do here is change the label for % the textboxes that ask for the initial conditions so it is clear to % the user what to enter. case 'pert' cancel = false; if (get(myfig.check.pert,'Value') == 1) set(myfig.label.Uic,'string','e(x,t)'); set(myfig.label.Sic,'string','d(x,t)'); else set(myfig.label.Uic,'string','u(x,t)'); set(myfig.label.Sic,'string','s(x,t)'); end return; % The radio buttons in the figure specify whether the user will define % the viewing window on the axes, or if Matlab will. If this is left % up to Matlab, then the textboxes where the user may enter the values % for the veiwing window are automatically disabled. Otherwise, % they're enabled. case 'axis' cancel = false; if (get(myfig.radio.top.auto,'Value') == 1) set(myfig.text.topymin,'enable','off'); set(myfig.text.topymax,'enable','off'); else set(myfig.text.topymin,'enable','on'); set(myfig.text.topymax,'enable','on'); end if (get(myfig.radio.bottom.auto,'Value') == 1) set(myfig.text.bottomymin,'enable','off'); set(myfig.text.bottomymax,'enable','off'); else set(myfig.text.bottomymin,'enable','on'); set(myfig.text.bottomymax,'enable','on'); end return; % When the user clicks on a specified object during the excecution of a % calculation. This sends the 'button' command to this function. % The calculation will stop by setting cancel = true. Appropriately % the status bar will display "Operation Cancelled". case 'button' status('Operation Cancelled'); progress(''); cancel = true; return; % One of the checkboxes in the figure corresponds to whether the % the viscous term will be considered in the PDE. When the checkbox % is clicked the 'visco' command is sent to this file. If the checkbox % is on then the textbox where the user may specify the value of % epsilon will be enabled, otherwise disabled. case 'visco' if (strcmp(get(myfig.text.eps,'enable'),'off')) set(myfig.text.eps,'enable','on','visible','on','string','0'); else set(myfig.text.eps,'enable','off','visible','off'); end return; case 'analyze' set(myfig.fig2,'visible','on'); if (get(myfig.figure2.menu,'Value') == 1) analyze(); else characteristics(); end return; case 'menu2' if (get(myfig.figure2.menu,'Value') == 2) set(myfig.figure2.measure,'string','0.25'); set(myfig.figure2.mlabel,'string','Refinement'); set(myfig.figure2.measure2,'Visible','off'); set(myfig.figure2.m2label,'Visible','off'); set(myfig.figure2.checkL,'Visible','on'); set(myfig.figure2.checkR,'Visible','on'); else set(myfig.figure2.mlabel,'string','Measure'); set(myfig.figure2.measure,'string','int((L-Lw)^2+(R-Rw)^2)'); set(myfig.figure2.measure2,'string','0'); set(myfig.figure2.measure2,'Visible','on'); set(myfig.figure2.m2label,'Visible','on'); set(myfig.figure2.checkL,'Visible','off'); set(myfig.figure2.checkR,'Visible','off'); end return; case 'examples' examples(get(myfig.menu3, 'Value')-1); return; end end % The settings() function is called when the user decided to run a brand % new simulation. All of the user specified information is retrieved in % this function and stored in the appropriate structures. function settings() global myfig; global grid; global parameters; global playopts; global cancel; status('Setting up solution grid...'); if (cancel == true) return; end % All of the important user specified variables pertaining to running % the simulation are retrieved here. % % xmin := the lower spatial boundary of the solution's domain. % xmax := the upper spatial boundary of the solution's domain. % tfinal := the upper bound of the time domain for the solution % a := the parameter alpha in the original PDE. % b := the parameter beta in the original PDE. % uL := the left asymptotic boundary condition for u(x,t). % uR := the right asymptotic boundary condition for u(x,t). % xprec := 1/dx, the number of meshpoints per x unit. xmin = str2num(get(myfig.text.xmin,'string')); xmax = str2num(get(myfig.text.xmax,'string')); tfinal = str2num(get(myfig.text.tfinal,'string')); a = str2num(get(myfig.text.alpha,'string')); b = str2num(get(myfig.text.beta,'string')); uL = str2num(get(myfig.text.uL,'string')); uR = str2num(get(myfig.text.uR,'string')); xprec = str2num(get(myfig.text.xprec,'string')); % Now we store the values pertaining to the solution mesh in the grid % structure. grid.xInterval = [xmin xmax]; grid.tInterval = [0 tfinal]; grid.xmin = xmin; grid.xmax = xmax; grid.tfinal = tfinal; grid.xprec = xprec; % In some places it will be easier to use dx rather than xprec, so this % is also stored. Also, the number of x meshpoints is calculated and % stored. grid.dx = 1/grid.xprec; grid.X = floor(grid.xprec*(xmax - xmin)); % NGC is the number of ghost cells defined in the mesh. This will be % necessary in our case since we would like the spatial derivative of u % on the boundaries to be 0. grid.NGC = 2; % We also store the partition of the grid as a vector in the grid % structure. grid.x = (xmin):(grid.dx):(xmax); % We'll also need to be able to easily distinguish between the actual % spatial domain (not including ghost cells), and the virtual spatial % domain (including ghost cells). The corresponding indices for the % upper and lower boundaries of these two intervals are calculated and % stored in the grid structure. grid.x0 = grid.NGC+1; grid.xN = grid.NGC+1+grid.X; grid.X0 = 1; grid.XN = grid.X+1+grid.NGC*2; % The essential PDE parameters are logically stored in the parameters % structure. parameters.a = a; parameters.b = b; parameters.uL = uL; parameters.uR = uR; parameters.pert = get(myfig.check.pert,'Value'); parameters.viscoterm = get(myfig.check.viscoterm,'Value'); % If the user indicates that they wish to include the viscous term in % the PDE, then the desired value of epsilon is retreived and stored % in the parameters structure. if (parameters.viscoterm) parameters.eps = str2num(get(myfig.text.eps,'string')); else parameters.eps = 0; end % In the playopts structure, there is an option for the frame rate. % This is set to be fixed at 10 frames/second. playopts.fps = 50; % This line tells the program to always calculate the wave solution, % regardless of whether or not the user wishes it to be displayed. parameters.wave = true; end % clearSolution() simply clears the four matrices of the sol structure. function clearSolution() global sol; global sol2; sol.u = []; sol.s = []; sol.uw = []; sol.sw = []; sol2.u = []; sol2.s = []; sol2.uw = []; sol2.sw = []; end % initialConditions() uses the user input to compute the initial condition % for the simulation on the already defined solution mesh. The function % returns a structure containing four columb vectors which is then intended % to be copied into the first row of the solution matries. The initial % condition is calculated for the function u(x,t), s(x,t), and the initial % condition for the corresponding traveling wave solution given the user % specified parameter set. function ic = initialConditions() global myfig; global grid; global sol; global sol2; global parameters; global cancel; % If the user cancels the process we return set the matrices to zero to % indicate the user cancelation rather than returning an set of empty % matrices which would indicate a calculation error on return. if (cancel == true) ic.u = [0]; ic.s = [0]; ic.uw = [0]; ic.sw = [0]; return; end % If we have already run a simulation but are simply continuing in time % on the same simulation, then all we need to do in this function is % take the last columb of the previous solution matrices and use those % as our new initial conditions. Then we're done so the function % returns. if (~isempty(sol.u)) [m n] = size(sol.u); ic.u = sol.u(:,n); ic.s = sol.s(:,n); if (parameters.wave) ic.uw = sol.uw(:,n); ic.sw = sol.sw(:,n); end return; end if (~isempty(sol2.u)) [m n] = size(sol2.u); ic.u = sol2.u(:,n); ic.s = sol2.s(:,n); if (parameters.wave) ic.uw = sol2.uw(:,n); ic.sw = sol2.sw(:,n); end return; end % If we're running a new simulation then we must proceed from here. % The following definitions retreive data from the global structures % grid and paramters, and simply store them in variables that are more % concisely named. x = grid.x; xInterval = grid.xInterval; X0 = grid.X0; x0 = grid.x0; xN = grid.xN; XN = grid.XN; NGC = grid.NGC; a = parameters.a; b = parameters.b; uL = parameters.uL; uR = parameters.uR; c = (uL+uR)/2; theta = uL*uR/2; pert = parameters.pert; % The status bar is updated to indicate at what point in the % calculation we have come to. status('Computing Initial Condition...'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % If the user has specified that they wish to view the traveling wave % solution as well, then we must also provide an initial condition for % the traveling wave solution (above this was set to always be true). % In fact this option is also needed if the user specifies the initial % condition as a wave pertubation instead of a standalone function. if (parameters.wave) [sw uw] = wavesolution(a,b,uL,uR,xInterval,x); end % The raw mathematical text that the user has entered to describe the % initial conditions are copied into the variables uic and sic. uic = get(myfig.text.Uic,'string'); sic = get(myfig.text.Sic,'string'); if (strcmp(sic,'WAVE') | strcmp(sic,'NEWT')) auto = true; else auto = false; end % The text is now parsed using the parse() function which is described % in greater detail below. The valuable data is then stored in the % init structure. The parse() function is necessary to discen % piecewise functions. In the text field the user may enter a function % in terms of x followed by $a,b$ which means the function is only % defined on the interval [a,b]. Several function may follow this % symbol followed by another interval on which that function is % defined. Outisde of the specified intervals, the function is assumed % to be zero. The parse() function stores this data in the init % structure in such a way that makes it easy to loop through to % calculate the actual function values in terms of x along the solution % domain. init = parse(uic,sic); % Here we initialize the two initial condition vectors for u(x,t) and % s(x,t). u = zeros(size(x)); s = zeros(size(x)); % Here we extract the necessary information from the init structure and % calculate the intended initial condition specified by the user. One % loop is needed for each function. This process is tried because in % case of an error we'd prefer it be displayed in the status bar and % not as the Matlab output. If an error does occur we set the initial % condition matrices to empty, as promised above. try if (init.u.size ~= 0) % This loop goes through each piecewise component of the the % function, evaluates it on the specified interval and adds it % to the current function that has been calculated in previous % steps. At the end of the loop we end up with the intended % piecewise function defined on the user specified grid. This % is stored in the vectors u and s respectively. for i = 1:(init.u.size) if (init.u.llim(i) == x(1)) u = u + eval(init.u.exp(i).str).*(x >= init.u.llim(i)& x <= init.u.rlim(i)); else u = u + eval(init.u.exp(i).str).*(x > init.u.llim(i) & x <= init.u.rlim(i)); end end else % In the event that the function was not piecewise, then our % job is much easier and no loop is necessary. u = eval(init.u.exp(1).str).*ones(size(x)); end % Finally we need to be sure that the dimensions of u match those % of x, where x is the vector containing the partition of the % solution grid. u = u.*ones(size(x)); % Repeat the analogous process for the function s(x,t) if (~auto) if (init.s.size ~= 0) for i = 1:(init.s.size) if (init.s.llim(i) == x(1)) s = s + eval(init.s.exp(i).str).*double(x >= init.s.llim(i) & x <= init.s.rlim(i)); else s = s + eval(init.s.exp(i).str).*double(x > init.s.llim(i) & x <= init.s.rlim(i)); end end else s = eval(init.s.exp(1).str).*ones(size(x)); end s = s.*ones(size(x)); end catch % If an error occurs the matricies of the ic structure are set to % empty so that the error is detected in the parent function. disp('We caught an error!!') ic.s = []; ic.u = []; ic.uw = []; ic.sw = []; return; end % If the user has selected the option of defining the initial condition % as a perturbation of the wave as opposed to a standalone function, % then we must add the wave solution to the calculated initial % condition, which is actually just a perturbation. if (pert & parameters.wave & ~auto) s = s + sw; u = u + uw; end if (pert & parameters.wave & auto) u = u + uw; if (strcmp(sic,'WAVE')) s = u.^2/2-c*u+theta; end if (strcmp(sic,'NEWT')) s = (a^2/b*dx(u'))'; end end area = int((u-uw)'); translation = area/(uL-uR); if (translation < 0) xInterval = [grid.xmin grid.xmax-translation]; else xInterval = [grid.xmin-translation grid.xmax]; end [sw uw] = wavesolution(a,b,uL,uR,xInterval,grid.x-translation); % In certain cases, we find that we cannot calculate a wave solution % with the parameters settings the user has defined. In this case, % parameters.wave would be set to false and hence we will return an % error saying that a perturbation cannot be defined because a wave % solution cannot be defined. if (pert & ~parameters.wave) errmsg('No wave solution exists for given parameter settings!'); return; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % After defining the initial conditions correctly for both u(x,t) and % s(x,t), we must then define the ghost cells correctly. As stated, we % assume the spatial derrivative is zero on the boundaries. Thus we % need the ghost cells on each end to match the values on either end of % the initial condition. So the vectors u and s are widended % appropriated and filled in with the correctly valued ghost cells. [m n] = size(x); u = [ones(1,NGC)*u(1) u ones(1,NGC)*u(n)]; s = [ones(1,NGC)*s(1) s ones(1,NGC)*s(n)]; % If the user wishes to view a wave solution, then the initial % condition vectors for this must be adjusted analogously. if (parameters.wave) uw = [ones(1,NGC)*uw(1) uw ones(1,NGC)*uw(n)]; sw = [ones(1,NGC)*sw(1) sw ones(1,NGC)*sw(n)]; end % We want to pass a comlumb vector, not a row vector, from this % function. That is because in the solution matrices, each row % represents a fixed point in space, and each columb represents a fixed % time, so the initial condition should be a comlumb vector. Thus the % two vectors are transposed and stored in the ic structure. ic.u = u'; ic.s = s'; % If the user wishes to view the wave solutions, then the analogous % thing is done for the initial condition vectors for the wave % solutions. If the user does not wish to see the wave solutions, then % these vectors are set to empty in order to conserve memory for the % program, to avoid unnecessary trouble later on. if (parameters.wave) ic.uw = uw'; ic.sw = sw'; else ic.uw = []; ic.sw = []; end % One may be curious to ask why the user is allowed to specify the % spatial grid precision but not the time grid precision. In order to % meet the necessary stability criterion, we need to make sure that dt % is small enough, or that the t-precision is large enough. When % introducing a viscous term into the system, dt needs to be on the % order of dx^2. On the other hand, if the viscous term is not % introduced, then t-precision must be greater than a certain constant % multiple of x-precision. This constant multiple is calculated and % stored as the variable ratio, then used to calculate the t-precision. % Finally, the new important data concerning the mesh with regard to % time is stored in the grid structure. if (parameters.viscoterm) grid.tprec = (grid.xprec)^2; else ratio = max(abs(ic.u))+parameters.a+1/parameters.b; grid.tprec = floor(grid.xprec*ratio)+1; end grid.dt = 1/grid.tprec; grid.T = floor(grid.tfinal*grid.tprec); end % simulateSolution() is pretty much the heart and soul of this program. It % does all of the actual PDE simulations and everything else about this % program is to help the appearence, either graphically or in the code. % The function takes in two arguments. The first is the structure ic which % contains information about the initial conditions that were calculated % using initialCodition(). The second argument finiteDiff is a boolean % telling which method should be used to simulate the solution of the PDE. % If true we use the Lax-Wendroff method, and if false we'll use a finite % volume scheme. % simulateSolution() takes up more memory and time than any other function % does since it's doing the hard calculations. It runs through a loop, % where each cycle corresponds to one point on the time grid, which can % multiply to become a very large number of cycles. % The function contiually updates the status bar and time bar so that the % user has a sense of how much longer the calculation will take. Also, it % checks after each cycle whether the user has cancelled the calculation. % After the loop is finished, the amount of memory being used is displayed % in the status bar so that the user may compare it to the amount of unused % memory on the system. Finally, the solution matrices, containing the % function values at each gridpoint in the mesh, are stored in the sol % structure and the sol structure is returned by the function. function solution = simulateSolution(ic,method) global parameters; global grid; global myfig; global playopts; global cancel; % Update the status bar. status('Running Simulation...'); if (cancel == true) return; end % Retrieve the parameters from the various global structures and store % them in variables having more concise names for easier readability % later on. The new variables, not already defined in settings(), are % defined below. % % T := the number of columbs in the solution matrix, or otherwise, % the number gridpoints corresponding to time. % % dt := the spacing between grid points along the time axis. % % tprec: the number of grid points per unit time, along the time % axis. % % eps := epsilon, the coeficient on the viscous term of the PDE. x = grid.x; xInterval = grid.xInterval; X0 = grid.X0; x0 = grid.x0; xN = grid.xN; XN = grid.XN; T = grid.T; dt = grid.dt; dx = grid.dx; NGC = grid.NGC; tprec = grid.tprec; xmin = grid.xmin; xmax = grid.xmax; a = parameters.a; b = parameters.b; uL = parameters.uL; uR = parameters.uR; eps = parameters.eps; fps = playopts.fps; % The initial conditions are extracted from the ic structure and each % vector is stored as two different variables. So two copies are made % of each. This is because of the way that the data will be compressed % as the function is being calculated. While it certainly increases % accuracy to consider time on a very fine grid, it certainly doesn't % help the appearence during playback that much more. So only the most % important instants in time are retained while the rest are deleted % from memory after being used in the simulation. This greatly reduces % the amount of memory used by this function and program in general, % the sol structure store more data than anything else. % One of the copies ends up storing all of the necessary snapshots in % time, while the other serves as a placeholder for the new initial % condition after the other time values are erased from memory. uu = ic.u; ss = ic.s; u = ic.u; s = ic.s; if (parameters.wave) uuw = ic.uw; ssw = ic.sw; uw = ic.uw; sw = ic.sw; end % In the event that the user wishes to use the finite volume method, it % will be easier to think of one vector q containing a component for u % and another for s, instead of thinking of u and s as two seperate % functions. if (method==2) q = [ic.u' ; ic.s']; if (parameters.wave) qw = [ic.uw' ; ic.sw']; end end if (method == 3 | method == 4) LLxi = ic.u' - 2*sqrt(a^2 + ic.s'); RReta = ic.u' + 2*sqrt(a^2 + ic.s'); LLeta = LLxi; RRxi = RReta; L = LLxi; R = RReta; bcL = [L(x0) L(xN)]; bcR = [R(x0) R(xN)]; lambdaXi = zeros(size(L)); lambdaEta = zeros(size(R)); xXi = x; xEta = x; for k=1:NGC xXi = [xmin-k*dx, xXi, xmax+k*dx]; xEta = [xmin-k*dx, xEta, xmax+k*dx]; end dLdxi = zeros(size(L)); dRdeta = zeros(size(R)); if (parameters.wave) LLxiw = ic.uw' - 2*sqrt(a^2 + ic.sw'); RRetaw = ic.uw' + 2*sqrt(a^2 + ic.sw'); LLetaw = LLxiw; RRxiw = RRetaw; Rw = RRetaw; Lw = LLxiw; xXiw = x; xEtaw = x; for k=1:NGC xXiw = [x(1)-k*dx, xXiw, x(xN-NGC)+k*dx]; xEtaw = [x(1)-k*dx, xEtaw, x(xN-NGC)+k*dx]; end end xFull = [ones(1:NGC)*xmin, x, ones(1:NGC)*xmax]; if (method == 4) lambdaLxi = zeros(size(L)); lambdaReta = zeros(size(R)); dLdxixi = zeros(size(L)); dRdetaeta = zeros(size(R)); end end % Next we define a few variables that contain information to be reused % again and again during the loop. % % i := the range of x indices of the solution grid not including the % ghost cells. % % z := a counter which recycles every 200 loops so that we know when to % update the time bar, as is done periodically. i = x0:xN; z=0; % There are several values that are used as coefficients over and over % again in each loop, so rather than having the computer take the time % to calculate them everytime, we simply calculate them once here and % store them in the following variables. Note that c4 is in fact a % vector, this is used to make sure that the solution vectors % calculated below have the correct dimensions. c1 = dt/dx; c2 = c1^2/2; c3 = exp(-b*dt); c4 = ones(NGC+1,1); % These set of vectors will be used to store finite difference % calculations for each point along the spatial grid. This needs to be % repeated at each time, but we've optimized the program so that it % only needs to carry these calculations out once per time step, rather % than several. u0 = zeros(size(i)); u1 = zeros(size(i)); u2 = zeros(size(i)); u12 = zeros(size(i)); s0 = zeros(size(i)); s1 = zeros(size(i)); s2 = zeros(size(i)); % Here are several more variables that will serve a few more purposes % during the execution of this loop. % % t1 := the current time marked before the loop begins. % % frame := the frame counter, refering to playback frames. % % t := the counter for the time grid point in the current frame. This % value is recycled when beginning to construct a new frame. % % memusage := a structure that is initialized to store information % about the memory being used by this calculation. t1=clock; frame = 1; t = 0; memusage = []; for n=1:T time = dt*n; t = t+1; % It is important now to check if the user wishes to cancel before % each loop, otherwise there will be a very large time delay % between when the user wishes to cancel, and when the program % actually cancels. if (cancel == true) return; end % Depending on whether finiteDiff is true or false, we may execute % one of two blocks of code on each loop. switch (method) case 1 % Use the first order finite differences to approximate the % values of the derivatives at each point for these functions. % All of the derrivative approximations are in terms of x and % pertain to the current time. % % u0(x) := u(x,t) % u1(x) := u'(x,t) % u2(x) := u''(x,t) % u12(x) := (u^2(x,t))' % s0(x) := s(x,t) % s1(x) := s'(x,t) % s2(x) := s''(x,t) u0(i) = u(i,t); u1(i) = (u(i+1,t) - u(i-1,t))/2; u2(i) = u(i+1,t) - 2*u(i,t) + u(i-1,t); u12(i) = (u(i+1,t).^2-u(i-1,t).^2)/4; s0(i) = s(i,t); s1(i) = (s(i+1,t) - s(i-1,t))/2; s2(i) = s(i+1,t) - 2*s(i,t) + s(i-1,t); % Now this is where the chosen algorithm itself is actually % executed. The next time step is calculated based on the % current timestep using the Lax-Wendroff method. The % functions ut, utt, st and stt approximate the time % derrivatives of these function in terms of the x derrivatives % of the corresponding functions. % The viscous relaxation is multiplied in at the end of the % calculation and the effect of the viscous term is added in % after that. u(i,t+1) = u0(i) + c1*ut(s1(i),u12(i)) + ... c2*utt(u0(i),u1(i),u2(i),s0(i),s1(i),s2(i),a); s(i,t+1) = s0(i) + c1*st(u0(i),u1(i),s0(i),s1(i),a) + ... c2*stt(u0(i),u2(i),s0(i),s1(i),s2(i),a); s(i,t+1) = c3*s(i,t+1); u(i,t+1) = u(i,t+1)+eps*u2(i)'; % Next we have to make sure that the ghost cells are filled in. % They need to be set to the same value as the values on the % extreme ends of the interval, so that the x derivative of the % boundaries is always 0. u(1:x0,t+1) = c4*u(x0,t+1); u(xN:XN,t+1) = c4*u(xN,t+1); s(1:x0,t+1) = c4*s(x0,t+1); s(xN:XN,t+1) = c4*s(xN,t+1); % If the user wishes to display the wave solution as well, then % then the wave solution is simulated analogously below. if (parameters.wave) u0(i) = uw(i,t); u1(i) = (uw(i+1,t) - uw(i-1,t))/2; u2(i) = uw(i+1,t) - 2*uw(i,t) + uw(i-1,t); u02(i) = (uw(i+1,t).^2-uw(i-1,t).^2)/4; s0(i) = sw(i,t); s1(i) = (sw(i+1,t) - sw(i-1,t))/2; s2(i) = sw(i+1,t) - 2*sw(i,t) + sw(i-1,t); uw(i,t+1) = u0(i) + c1*ut(s1(i),u02(i)) + ... c2*utt(u0(i),u1(i),u2(i),s0(i),s1(i),s2(i),a); sw(i,t+1) = s0(i) + c1*st(u0(i),u1(i),s0(i),s1(i),a) + ... c2*stt(u0(i),u2(i),s0(i),s1(i),s2(i),a); sw(i,t+1) = c3*sw(i,t+1); uw(i,t+1) = uw(i,t+1)+eps*u2(i)'; uw(1:x0,t+1) = c4*uw(x0,t+1); uw(xN:XN,t+1) = c4*uw(xN,t+1); sw(1:x0,t+1) = c4*sw(x0,t+1); sw(xN:XN,t+1) = c4*sw(xN,t+1); end % The below block of code will be executed at each cycle of this % loop if finiteDiff is false. The user has indicated that they'd % rather use a finite volume method to simulate the solution of % this PDE. case 2 % Begin by filling in the ghost cells correctly. for j=1:NGC q(:,x0-j) = q(:,x0); q(:,xN+j) = q(:,xN); end % Then run the simulation accounting for advection. The % advect() function will be described in much greater detail % below. q=advect(q,dt); % Now that q has been modified it is necessary to update the % ghost cells once again. for j=1:NGC q(:,x0-j) = q(:,x0); q(:,xN+j) = q(:,xN); end % Now run the simulation accounting for the elastic effect in % the fluid. elastic() function will be described in much % greater detail below. q=elastic(q,dt); % Account for the relaxation of stress in the fluid % by multiplying the values of s by an exponentially decaying % constant. q(2,:) = q(2,:).*c3; % Finally, store the calculated solution in the slot for the % next time step in the matrices u and s respectively. u(:,t+1) = q(1,:)'; s(:,t+1) = q(2,:)'; % If the user wishes to view the wave solution as well, then % simulate the traveling wave solution analogously. if (parameters.wave) for j=1:NGC qw(:,x0-j) = qw(:,x0); qw(:,xN+j) = qw(:,xN); end qw=advect(qw,dt); for j=1:NGC qw(:,x0-j) = qw(:,x0); qw(:,xN+j) = qw(:,xN); end qw=elastic(qw,dt); qw(2,:) = qw(2,:).*c3; uw(:,t+1) = qw(1,:)'; sw(:,t+1) = qw(2,:)'; end % This is a first order numerical method of characteristics. % We have solved the Riemann Invariants of the original PDE and % are now using the characteristic information to numerical % solve the transformed PDE. At each time step the wave speed % is first calculated so that the characteristics can be % updated appropriately over one unit of time. Then, using the % source term the transformed system, we can calculate how L % and R change along each of these characteristics. After % calculating the updated L and R values, we interpolate the % solution back onto the original grid to eliminate having to % have the two characteristics meet at a particular point. case 3 % A few of these will need to be initialized above lambdaXi(i) = (3*RRxi(i)+LLxi(i))/4; lambdaEta(i) = (RReta(i)+3*LLeta(i))/4; dLdxi(i) = b*(RRxi(i) - LLxi(i))/4 - 4*a^2*b./(RRxi(i) - LLxi(i)); dRdeta(i) = -b*(RReta(i) - LLeta(i))/4 + 4*a^2*b./(RReta(i) - LLeta(i)); % Update the sample points on the x-grid for the two % characteristics xXi(i) = xXi(i)+lambdaXi(i)*dt; xEta(i) = xEta(i)+lambdaEta(i)*dt; % The value of L and R along the xi and eta characteristics % respectively LLxi(i) = LLxi(i) + dLdxi(i)*dt; RReta(i) = RReta(i) + dRdeta(i)*dt; % To handel the boundary cases, we choose points far out of the % range of influence. We map where the characteristics are at % and choose points that are at least three times the width of % the range to the left or right of the endpoints. Because L % and R are invariant along both characteristics at their % boundary values, then we can assume these points are equal to % the boundary values. farLeft = xEta(x0); farRight = xXi(xN); range = farRight-farLeft; for k=1:NGC xXi(x0-k) = farLeft - 3*range - k*dx; xEta(x0-k) = farLeft - 3*range - k*dx; xXi(xN+k) = farRight + 3*range + k*dx; xEta(xN+k) = farRight + 3*range + k*dx; end % Now that the boundaries have been accounted for, we can % interpolate properly along the opposite characteristics. LLeta(i) = interp1(xXi,LLxi,xEta(i)); RRxi(i) = interp1(xEta,RReta,xXi(i)); % If the user wishes to display the wave solution as well, then % then the wave solution is simulated analogously below. if (parameters.wave) lambdaXi(i) = (3*RRxiw(i)+LLxiw(i))/4; lambdaEta(i) = (RRetaw(i)+3*LLetaw(i))/4; dLdxi(i) = b*(RRxiw(i) - LLxiw(i))/4 - 4*a^2*b./(RRxiw(i) - LLxiw(i)); dRdeta(i) = -b*(RRetaw(i) - LLetaw(i))/4 + 4*a^2*b./(RRetaw(i) - LLetaw(i)); xXiw(i) = xXiw(i)+lambdaXi(i)*dt; xEtaw(i) = xEtaw(i)+lambdaEta(i)*dt; LLxiw(i) = LLxiw(i) + dLdxi(i)*dt; RRetaw(i) = RRetaw(i) + dRdeta(i)*dt; farLeft = xEta(x0); farRight = xXi(xN); range = farRight-farLeft; for k=1:NGC xXiw(x0-k) = farLeft - 3*range - k*dx; xEtaw(x0-k) = farLeft - 3*range - k*dx; xXiw(xN+k) = farRight + 3*range + k*dx; xEtaw(xN+k) = farRight + 3*range + k*dx; end LLetaw(i) = interp1(xXiw,LLxiw,xEta(i)); RRxiw(i) = interp1(xEtaw,RRetaw,xXiw(i)); end % A second order numerical method of characteristics shcheme. % This incorperates the Lax-Wendroff method for calculating the % movement of the charateristic and the change in L and R up to % second order accuracy. end % After the solution has been calculated for the next time step, it % may be necessary to compress memory depending on precisely which % time step we're on. We only want to keep time steps that % correspond to a frame snapshot during playback. Since we've % decided our playback speed to be 10 frames/sec, then we'll only % need to keep 10 evenly spaced time steps. % If this so happens to be one of the specially labeled time steps % then the following block of code is executed in order to conserve % memory. We take the calculation for the current time step, and % store it as the next frame in one of the vectors uu, ss, uuw, or % ssw. Then we erase the matrices u, s, uw and sw but then % initialize their first columbs as the frame just added. Then % they act as the new initial condition above. We increment the % frame counter by one, and recycle the number of time steps that % have occured since the saving of the last frame. if (n == floor(frame*tprec/fps)) frame = frame+1; if (method ~= 3 & method ~= 4) [trash sz] = size(u); uu(:,frame) = u(:,sz); ss(:,frame) = s(:,sz); uuw(:,frame) = uw(:,sz); ssw(:,frame) = sw(:,sz); u = uu(:,frame); s = ss(:,frame); uw = uuw(:,frame); sw = ssw(:,frame); end if (method == 3) Lrow = interp1(xXi,LLxi,xFull); Rrow = interp1(xEta,RReta,xFull); Lwrow = interp1(xXiw,LLxiw,xFull); Rwrow = interp1(xEtaw,RRetaw,xFull); L = [L ; Lrow]; R = [R ; Rrow]; Lw = [Lw ; Lwrow]; Rw = [Rw ; Rwrow]; end t=0; end % At the 100th cycle of the calculation, the current time is marked % and compred to the time from before the loop began. This % difference in time is then divided by 100 to get the average time % per cycle and then multiplied by the total number of cycles that % are expected to run. By this formula we have estimated the time % that the entire calculation will take. We then display this time % in the status bar. if (n==100) t2 = clock; if (t1(5) ~= t2(5)) time = 60-t1(6)+t2(6); else time = t2(6)-t1(6); end esttime = T*time/100; msg = ['Running Simulation...',' Estimated Time =',num2str(esttime), ' seconds']; status(msg); end % Every 200 cycles we will update the time bar so that it displays % the current percentage of progress that has been made based on % the number of cycles that have occured. We know that the index n % must count from 1 through T. So the percentage n/T*100 is what % marks the progress that has been made. This percentage is % displayed in the time bar. z=z+1; if (z == 200) progress = n/T*100; msg = ['Simulation is ',num2str(progress,2),'% complete.']; set(myfig.progress,'string',msg); getframe(gcf); z=0; end % mem = whos; % [numb trash] = size(mem); % memory = 0; % for q = 1:numb % memory = memory + mem(q).bytes; % end % memusage = [memusage ; memory]; %disp(['Using ',num2str(memory/1000000,3),'MB of memory.']); end % maxmemusage = max(memusage); % disp(['Max memory usage = ',num2str(maxmemusage/1000000,3),'MB']); % The solution that were calculated in the above loop are now stored in % sol structure, which is the output of this function. if (method == 3 | method == 4) uu = (L+R)/2; ss = (L-R).^2/16-a^2; uuw = (Lw+Rw)/2; ssw = (Lw-Rw).^2/16-a^2; uu = uu'; ss = ss'; uuw = uuw'; ssw = ssw'; end solution.u = uu; solution.s = ss; if (parameters.wave) solution.uw = uuw; solution.sw = ssw; end % Then the progress bar is updated to indicate that the simulation is % complete. set(myfig.progress,'string','Simulation complete.'); getframe(gcf); end % play() is the function that is called when we wish to play back the % solution that has been calculated. Presumably the solution has been % calculated out from time=0 until time=T where T is some user specified % ending time. The user may not wish to playback the entire length of the % calculated solution, but instead only a subset of the total time segment. % Thus, the two input arguments for play indicate the start and finish % time of that segment. % The function first retreives all necessary information from the global % structures. It converts the start/finish times to points on the time % grid. Then it checks how the user has defined the axes and sets up to % plot the solution according to these definitions. Then the user % specified viewing window for the axes are set up, if the user chose to % specify them manually. Then we run a loop which sequentially displays % the condensed calculated solution at the desired frame rate. Once this % process is complete the total memory usage by the program is calculated % and displayed in the status bar. function play(start,finish) global playopts; global myfig; global grid; global parameters; global sol; global sol2; global cancel; % Indicate to the user that the solution will begin playing % momentarily. status('Playing Solution...'); if (cancel == true) return; end playopts.start = start; % Here we extract the necessary information contained in each structure % and store them as variables that are more concisely named. Most all % of the variables have been defined somewhere above. % % c := the speed of the traveling wave given the boundary conditions. x = grid.x; x0 = grid.x0; xN = grid.xN; NGC = grid.NGC; tprec = grid.tprec; xmin = grid.xmin; xmax = grid.xmax; T = grid.T; c = (parameters.uL + parameters.uR)/2; theta = parameters.uL*parameters.uR/2; a = parameters.a; b = parameters.b; fps = playopts.fps; save = get(myfig.check.save,'Value'); method1 = get(myfig.menu,'Value'); method2 = get(myfig.menu2,'Value')-1; meth2 = (method2 ~= 0 & method2 ~= method1); % A new structure called ax is created which contains the user % definitions of what will be displayed on each axes. Essentially the % user has the power to graph any function against any other. We store % these as strings because it is presumed the user will enter a % formula. ax.top.x = get(myfig.text.topx,'string'); ax.top.y = get(myfig.text.topy,'string'); ax.bottom.x = get(myfig.text.bottomx,'string'); ax.bottom.y = get(myfig.text.bottomy,'string'); % dispwave := true if the user wishes to view the traveling wave % solution along with the other solution. false if the user does not % wish this. dispwave = get(myfig.check.wave,'Value'); % In the event that a wave solution does not exist for the parameter % settings, it will not be possible to display the traveling wave % solution as the user would want, so an error message is displayed, % but we don't return because we can still playback the regular % solution. if (dispwave & ~parameters.wave) errmsg('A wave solution does not exist for the parameters settings!') end % Clearly we need 0 < start < finish. And of course finish must be % within the calculated time interval. if (start < 0 | start >= finish | finish > grid.tfinal) errmsg('Bad play time settings!'); return; end % Here we convert the start and finish time to indices corresponding to % the mesh we have defined. The function time2index() is discribed in % greater detail below. begin = time2index(start); final = time2index(finish); % If for some reason the finishing time was too large, we'll assume the % user just wishes to watch the movie until its end, so we'll simply % set the final index to the last on in the grid. if (final > T) final = T; end % i := a vector containing indices corresponding to the spatial % partition of the grid, not including ghost cells. % n := a vector containing indices corresponding to the time partition % of the grid within the desired time interval. i = x0:xN; n = begin:final; % x is the spatial partition itself, not including ghost cells. Since % the vector i will not start at 1, we need to offset it by the number % of ghost cells. This next step makes sure that the counter i matches % with the corresponding values in the paritition x. x(i) = x(i-NGC); % A solution of all zeros will be needed in case the axes definitions % the user has chosen leads to a trivial graph. sol.z = zeros(size(sol.u)); sol2.z = zeros(size(sol2.u)); % stringPrepAx() is defined in greater detail below but essentially it % converts the users inputted formula to matlab code that can be % evaluated. There are several shorthand symbols that are used in % place of code to make the user's life easier. % % u_ = u(x,t) % e_ = epsilon(x,t) = u(x,t) - uw(x,t) % s_ = s(x,t) % d_ = delta(x,t) = s(x,t) - sw(x,t) % % The derivatives of each of these functions may be calculated as well % using the dx(), dxx(), dt() and dtt() functions respectively. topx = stringPrepAx(ax.top.x,false); topy = stringPrepAx(ax.top.y,false); bottomx = stringPrepAx(ax.bottom.x,false); bottomy = stringPrepAx(ax.bottom.y,false); if (meth2) topx2 = stringPrepAx2(ax.top.x,false); topy2 = stringPrepAx2(ax.top.y,false); bottomx2 = stringPrepAx2(ax.bottom.x,false); bottomy2 = stringPrepAx2(ax.bottom.y,false); end % Now we calculate the actual values for each axes using the user % specified formulas. These take the solution matrices and apply the % user's formula to them. In the even that the eval() function cannot % evaluate the user's input, then clearly they have made a mistake so % an error is displayed in the status bar and the user must enter valid % axes definitions. try x1 = eval(topx); x2 = eval(bottomx); y1 = eval(topy); y2 = eval(bottomy); if (meth2) x12 = eval(topx2); x22 = eval(bottomx2); y12 = eval(topy2); y22 = eval(bottomy2); end catch errmsg('Bad axis setttings!'); return; end % Below is the analogous process to what was done above, only it is % done for the wave solution, but only if the user has chosen to have % the wave displayed. The second argument of stringPrepAx() simply % indicates whether the axes are being calculated for the wave solution % or not. if (dispwave & parameters.wave) topxw = stringPrepAx(ax.top.x,true); topyw = stringPrepAx(ax.top.y,true); bottomxw = stringPrepAx(ax.bottom.x,true); bottomyw = stringPrepAx(ax.bottom.y,true); if (meth2) topxw2 = stringPrepAx2(ax.top.x,true); topyw2 = stringPrepAx2(ax.top.y,true); bottomxw2 = stringPrepAx2(ax.bottom.x,true); bottomyw2 = stringPrepAx2(ax.bottom.y,true); end try xw1 = eval(topxw); xw2 = eval(bottomxw); yw1 = eval(topyw); yw2 = eval(bottomyw); if (meth2) xw12 = eval(topxw2); xw22 = eval(bottomxw2); yw12 = eval(topyw2); yw22 = eval(bottomyw2); end catch errmsg('Bad axis settings!!'); return; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [n1 trash] = size(x1); [n2 trash] = size(x2); clear trash; % Here we determine the smallest veiwing window needed to display the % solution on the set of axes for the time interval. In most cases % we're maximizing or minimizing a matrix, so the max/min function must % be called twice in this case to get the max/min of the entire matrix. if (n1 == 1) x1max = max(x1); x1min = min(x1); else x1max = max(max(x1(i,n))); x1min = min(min(x1(i,n))); end if (n2 == 1) x2max = max(x2); x2min = min(x2); else x2max = max(max(x2(i,n))); x2min = min(min(x2(i,n))); end % The analogous process is repeated for the wave solution, but only if % the user wishes the wave solution to be displayed. if (dispwave & parameters.wave) if (n1 == 1) xw1max = max(xw1); xw1min = min(xw1); else xw1max = max(max(xw1(i,n))); xw1min = min(min(xw1(i,n))); end if (n2 == 1) xw2max = max(xw2); xw2min = min(xw2); else xw2max = max(max(xw2(i,n))); xw2min = min(min(xw2(i,n))); end x1max = max(x1max,xw1max); x2max = max(x2max,xw2max); x1min = min(x1min,xw1min); x2min = min(x2min,xw2min); end % It may be the case that the max/min values will equal eachother % making the veiwing window infinitely small. In this case we widen it % by 0.1 units on either side. There's nothing special about the % number we have chosen here. if (x1min == x1max) x1max = x1min + 0.1; x1min = x1min - 0.1; end if (x2min == x2max) x2max = x2min + 0.1; x2min = x2min - 0.1; end % If the user has specified their own viewing window, then we obtain % the information for that viewing window here, otherwise we retain the % window we calculated above. if (get(myfig.radio.top.fixed,'Value') == 1) y1min = str2num(get(myfig.text.topymin,'string')); y1max = str2num(get(myfig.text.topymax,'string')); if (isempty(y1min) | isempty(y1max) | y1min >= y1max) errmsg('Invalid y-limit settings'); return; end else y1max = max(max(y1(i,n))); y1min = min(min(y1(i,n))); if (dispwave & parameters.wave) yw1max = max(max(yw1(i,n))); yw1min = min(min(yw1(i,n))); y1max = max(y1max,yw1max); y1min = min(y1min,yw1min); end if (y1min == y1max) y1max = y1min + 0.1; y1min = y1min - 0.1; end end % The analogous calculation is carried out for the second axes. if (get(myfig.radio.bottom.fixed,'Value') == 1) y2min = str2num(get(myfig.text.bottomymin,'string')); y2max = str2num(get(myfig.text.bottomymax,'string')); if (isempty(y2min) | isempty(y2max) | y2min > y2max) errmsg('Invalid y-limit settings'); return; end else y2max = max(max(y2(i,n))); y2min = min(min(y2(i,n))); if (dispwave & parameters.wave) yw2max = max(max(yw2(i,n))); yw2min = min(min(yw2(i,n))); y2max = max(y2max,yw2max); y2min = min(y2min,yw2min); end if (y2min == y2max) y2max = y2min + 0.1; y2min = y2min - 0.1; end end % Here we begin our main loop of the play cycle. The counter k, begins % at the first index of our time interval and ends at the last. We % take care of one set of axes at a time. The solution matrices, which % have now been acted on by the user specified formula to produce x1, % x2, y1, and y2 are used to plot the solution. At each cycle we plot % the kth row of x1 against that of y1 on the first axes. We do the % analogous thing for the second axes. Also, if the user wishes to see % the wave solution then these are plotted analogously as well. The % axes limits are reset (so that they don't change) after each cycle % and the graph is relabeled after each cycle. A short pause is taken % at the end of each cycle so the user may view the process as if it % were a movie. frame = 1; for k = begin:final % We need to check at the beginning of each cycle for a user % cancellation since this loop may run for a while in some cases. if (cancel == true) return; end % Here we take the current index and figure out which time it % corresponds to in the interval [0 tfinal]. The inxex2time() % function works in the opposite way as time2index() and is % described in greater detail below. time = index2time(k); % The first (upper) subplot is taken care of first. We make it the % current axes and then clear what was previously being displayed on % it. We turn hold ON in case we need to draw more than one thing % on the axes during this cycle. axes(myfig.Uax); cla; hold on; % Plot the kth columb of x1 against the kth columb of y1. If x1 is % defined simply to be the vector x, then it has only one columb % (i.e. n1 = 1) so there is no need to search for the kth columb. if (n1 == 1) plot(x1(i),y1(i,k),'-b'); if (meth2) plot(x12(i),y12(i,k),':b'); end else plot(x1(i,k),y1(i,k),'-b'); if (meth2) plot(x12(i,k),y12(i,k),':b'); end end % If the user wishes to display the wave solution as well then we'll % graph the corresponding columbs of xw1 against yw1. if (dispwave & parameters.wave) if (n1 == 1) plot(xw1(i),yw1(i,k),'-r'); if (meth2) plot(xw12(i),yw12(i,k),':r'); end else plot(xw1(i,k),yw1(i,k),'-r'); if (meth2) plot(xw12(i,k),yw12(i,k),':r'); end end end % Here we label the axes the same as the user specified formula, % only we erase every underscore. So instead of having u_ + s_ it % will simply display u + s. In the title we display the current % time, and of course this changes with every cycle so it creates % the effect of a timer as the solution is playing. Finally, the % axes limits are reset after each cycle so that they don't % accidentally change by Matlab!! This was the only way I found to % fix the axes limits throughout the loop. xlabel(strrep(ax.top.x,'_','')); ylabel(strrep(ax.top.y,'_','')); title(['| Time = ',num2str(time,2),' |']); axis([x1min x1max y1min y1max]); % The analogous steps are taken in constructing the second (lower) % subplot. axes(myfig.Sax); cla; hold on; if (n2 == 1) plot(x2(i),y2(i,k),'-b'); if (meth2) plot(x22(i),y22(i,k),':b'); end else plot(x2(i,n),y2(i,k),'-b'); if (meth2) plot(x22(i,n),y22(i,k),':b'); end end if (dispwave & parameters.wave) if (n2 == 1) plot(xw2(i),yw2(i,k),'-r'); if (meth2) plot(xw22(i),yw22(i,k),':r'); end else plot(xw2(i,n),yw2(i,k),'-r'); if (meth2) plot(xw22(i,n),yw22(i,k),':r'); end end end xlabel(strrep(ax.bottom.x,'_','')); ylabel(strrep(ax.bottom.y,'_','')); title(['| Time = ',num2str(time,2),' |']); axis([x2min x2max y2min y2max]); % A pause is taken so that the user may see the result of this cycle % before beginning a new one. This will have the effect of a flip % book, or cartoon. Thus the solution will appear animated in time % as it should. pause(0.001); if (save) Mov(frame) = getframe; frame = frame + 1; end end % Here we use the "whos" command to obtain information about the memory % currently being used by Matlab. We add the memory being used by each % variable together and display the total, in terms of Megabytes, in % the status bar. if (save) movie2avi(Mov, 'visco01.avi', 'compression', 'Indeo3', 'fps', 10); end mem = whos; [n trash] = size(mem); memory = 0; for i = 1:n memory = memory + mem(i).bytes; end status(['Using ',num2str(memory/1000000,3),'MB of memory.']); end % Before any simulations have been run, the user does not have the option % playing a solution back or continuing a simulation, because there is % nothing to playback or continue. Thus, these options are initially % disabled. Once the first simulation is run though, these options become % enabled to the user. The new options include the play and continue % buttons, the text field indicating the amount of time to continue % simulating, and the text fields indicating the start and stop time for % the playback. function activateButtons() global myfig; global grid; set(myfig.button.play,'enable','on'); set(myfig.text.playstart,'enable','on','string','0'); set(myfig.text.playend,'enable','on','string',num2str(grid.tfinal)); set(myfig.button.continue,'enable','on'); set(myfig.text.continue,'enable','on','string','5'); set(myfig.button.analyze,'enable','on'); end % The analyze() function is kept here JUST in case I might need it at some % point in the future, but is not essential to the actual program so I % won't go into great detail about it. Essentially it just opens up % another window and set of axes which displays whatever I have told it to % in the code below. It might plot a function measure against time for % example. This plot is not animated. function analyze() global sol; global parameters; global myfig; global grid; uL = parameters.uL; uR = parameters.uR; a = parameters.a; b = parameters.b; c = (uL + uR)/2; theta = uL*uR/2; measurestr = get(myfig.figure2.measure,'String'); measurestr = stringPrep(measurestr); measurestr = stringPrepAx(measurestr, false); y = eval(measurestr); measurestr2 = get(myfig.figure2.measure2,'String'); if (~strcmp(measurestr2,'')) str2 = true; measurestr2 = stringPrep(measurestr2); measurestr2 = stringPrepAx(measurestr2, false); y2 = eval(measurestr2); else str2 = false; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [m n] = size(y); t = index2time(1:n); % The other figure needs to be made visible at this time. % Set the figure's axes to the current set of axes. axes(myfig.figure2.ax); cla; hold on; % Plot both the measure over time and its upper bound which is a power % fucntion of time. plot(t,y,'-b'); if (str2) plot(t,y2,'-r'); end % Set the axes to fit around the measure's function over time. % if (max(y) <= 0) % axis([0 max(t) -0.1 0.1]); % else % axis([0 max(t) 0 max(y)]); % end axis tight; % Label the axes accordingly. xlabel('t'); ylabel('Energy'); end % The characteristics() function graphs the characteristics of the % solutions. function characteristics() global myfig; global sol; global parameters; global grid; refinement = str2num(get(myfig.figure2.measure,'String')); space = floor(refinement/grid.dx); u = sol.u; s = sol.s; a = parameters.a; b = parameters.b; dt = grid.dt; T = grid.tfinal; X = grid.X; dispLeft = get(myfig.figure2.checkL,'Value'); dispRight = get(myfig.figure2.checkR,'Value'); i = 1:space:X; x0 = grid.x(i); [m n] = size(x0); lambdaL = u + sqrt(s+a^2); lambdaR = u - sqrt(s+a^2); axes(myfig.figure2.ax); cla; hold on; xlabel('x'); ylabel('t'); options = []; for k = 1:n initdata = [x0(k)]; if (dispLeft) [tL,xL]=ode45('odeL',[0 T],initdata,options,lambdaL); plot(xL,tL,'-b'); end if (dispRight) [tR,xR]=ode45('odeL',[0 T],initdata,options,lambdaR); plot(xR,tR,'-r'); end axis tight; pause(0.001); end end % The advect() function was written by Professor Bob Guy and it basically % updates the current solution by accounting for the effect of advection. % The underlying philosophy behind Professor Guy's code is that the given % PDE can be decomposed into three smaller PDE's. One of them being a % purely advective PDE system shown below. % % function qnew=advect(q,dt); global grid; x0 = grid.x0; xN = grid.xN; dx = grid.dx; qnew = zeros(size(q)); % u is the solution to the Riemann problem % uedge = zeros(size(q(1,:))); % compute jumps for later use % du(x0-1:xN+2) = q(1,x0-1:xN+2) - q(1,x0-2:xN+1); ds(x0-1:xN+2) = q(2,x0-1:xN+2) - q(2,x0-2:xN+1); % solve Riemann problem % UL = q(1,x0-1:xN); UR = q(1,x0:xN+1); uedge(x0:xN+1) = double( (UL>0)&((UL+UR)>0) ).*UL + ... double( (UR<0)&((UL+UR)<0) ).*UR; % flux difference to update u % qnew(1,x0:xN) = q(1,x0:xN) ... - 0.5*dt/dx*( uedge(x0+1:xN+1).^2 - uedge(x0:xN).^2 ); % use edge velocities to propagate jumps is s % ds(x0:xN+1) = q(2,x0:xN+1) - q(2,x0-1:xN); qnew(2,x0:xN) = q(2,x0:xN) - ... dt/dx*( ds(x0:xN) .*uedge(x0:xN) .*double(uedge(x0:xN) >0)+ ... ds(x0+1:xN+1).*uedge(x0+1:xN+1).*double(uedge(x0+1:xN +1)<0)); end % function qnew=elastic(q,dt) global grid; global parameters; dx = grid.dx; x0 = grid.x0; xN = grid.xN; a = parameters.a; qnew = zeros(size(q)); % average s to the edges % sedge=zeros(size(q(2,:))); sedge(x0-1:xN+2) = 0.5 * ( q(2,x0-2:xN+1) + q(2,x0-1:xN+2)); % compute wave speeds at these edges % cp = sqrt(abs(sedge + a^2)); cm = -cp; % make right and left eigenvectors % Rp = [ones(size(cm)); cm]; Rm = [ones(size(cp)); cp]; Lp = [cm; ones(size(cm))]; Lm = [cp; ones(size(cp))]; [m i] = min(abs(dot(Lp,Rp))); % normalize left eigenvectors so that have biorthonormal set % Lp = Lp ./ [dot(Lp,Rp); dot(Lp,Rp)]; Lm = Lm ./ [dot(Lm,Rm); dot(Lm,Rm)]; % compute the jumps at the edges % dq(:,x0-1:xN+2) = q(:,x0-1:xN+2)-q(:,x0-2:xN+1); % compute strength of right and left going waves % ap = dot(Lp,dq,1); am = dot(Lm,dq,1); % for j=x0:xN % % qnew(:,j) = q(:,j) - cp(j)*dt/dx * Rp(:,j) * ap(j); % % qnew(:,j) = qnew(:,j) - cm(j+1)*dt/dx * Rm(:,j+1)* am(j+1); % % end j = x0:xN; qnew(1,j) = q(1,j) - (dt/dx)*cp(j).*Rp(1,j).*ap(j) - ... (dt/dx)*cm(j+1).*Rm(1,j+1).*am(j+1); qnew(2,j) = q(2,j) - (dt/dx)*cp(j).* Rp(2,j).*ap(j) - ... (dt/dx)*cm(j+1).*Rm(2,j+1).*am(j+1); F = zeros(size(q)); j = x0:xN+1; F(1,j) = F(1,j) + 0.5*abs(cp(j)).*( 1-dt/dx*abs(cp(j)) ).*Rp(1,j).*ap(j)... + 0.5*abs(cp(j)).*( 1-dt/dx*abs(cp(j)) ).*Rm(1,j).*am(j); F(2,j) = F(2,j) + 0.5*abs(cp(j)).*( 1-dt/dx*abs(cp(j)) ).*Rp(2,j).*ap(j) ... + 0.5*abs(cp(j)).*( 1-dt/dx*abs(cp(j)) ).*Rm(2,j).*am(j); j = x0:xN; qnew(1,j) = qnew(1,j) - dt/dx*( F(1,j+1)-F(1,j) ); qnew(2,j) = qnew(2,j) - dt/dx*( F(2,j+1)-F(2,j) ); end % This function uses the parameter settings to output the unique traveling % wave solution we found in our analysis. The solution is calculated on % the entire x-interval and the output is its interpolation on the % partition x. Internally it uses the functions U(p) and S(p) where % p = x - ct so that U(p) = u(x,t) and S(p) = u(x,t). Using this % substitution converts the system into an ODE which can be analyzed much % more simply. Since, in particular, we use this to solve for the initial % condition of the system, then really U(x) = u(x,0) and S(x) = s(x,0). function [sinit uinit] = wavesolution(a,b,uL,uR,xInterval,x) global xprec; global tprec; wave = 'weak'; I = max(abs(xInterval(1)),abs(xInterval(2))); P = a^2-(uL-uR)^2/4; d = abs(uL-uR); c = (uL+uR)/2; % Here we define the initial condition that will be plugged into ode45 to % solve for the wave equation. This point is the minimum S value and the % median value of U. U = c; S = -d^2/8; initdata = [U,S]; options = []; if (P >= -d^2/8 & P < 0) I = (d*sqrt(d^2+8*P)+8*P*atanh(sqrt(d^2+8*P)/d))/(2*b*d); end % if (P >= -d^2/8 & P < 0) % m = -d^2*b/(d^2+8*P) % I = abs(d/(2*m)) % yy = uL*(x<-I) + (c+m*x).*(x>=-I & x <= I) + uR*(x>I); % end if (P==0) I = d/(2*b); end if (P < -d^2/8) sinit = 0*(x ~= 0) - d^2/8*(x == 0); uinit = uL*(x <= 0) + uR*(x > 0); return; end % if (P >= -d^2/8 & P < 0) % m = -b*d^2/(d^2+8*P); % I = d/(2*m); % sinit % end % Starting from the initial point specified above, the solution runs both % fowards and backwards for the specified amount of time, given by the % xInterval. Backwards time is acheived simply by imputting -b instead of % b, since b is proportional to the time-scale. The system of equations % being solved for is described in the external function 'myode'. [tf,yf]=ode45('myode',[0 I],initdata,options,a,b,uL,uR); [tb,yb]=ode45('myode',[0 I],initdata,options,a,-b,uL,uR); % The backward solution needs to be reversed in order. Further, only one % value at 0 can be specified so this value is removed from the backward % solution. M = size(tb); M = M(1); tb = tb(2:M); yb = yb(2:M,:); M = M-1; % A is the permutation matrix which reverses the order of the vectors % defining the backward solution. A = zeros(M); for j = 1:M A(j,M-j+1)=1; end tb = -A*tb; yb = A*yb; % Combine the foward and backward soltutions into one large vector. t = [tb ; tf]; y = [yb ; yf]; % Interpolate the traveling wave solution on the partition defined by the % vector x. Output this interpolated solution which will be used as the % initial condition for this program's PDE algorithm. [m n] = size(x); sinit = zeros(size(x)); uinit = zeros(size(x)); for i = 1:n if (x(i)>-I & x(i)= I) sinit(i) = 0; uinit(i) = uL*(x(i) <= -I) + uR*(x(i) >= I); end end end % examples() % This function takes on a value specified by the examples menu selected by % the user. Depending on the selection this function modifies the input % fields accordingly so that the correct example is displayed when the user % runs the similation. The user still has the option of editing the input % values or selecting another example after this function is run. function examples(choice) global myfig; switch(choice) case 0 return; case 1 set(myfig.check.wave,'Value',1); set(myfig.check.pert,'Value',1); set(myfig.check.save,'Value',0); set(myfig.check.viscoterm,'Value',0); set(myfig.label.Uic,'string','e(x,t)'); set(myfig.label.Sic,'string','d(x,t)'); set(myfig.text.alpha, 'string','1') set(myfig.text.beta, 'string','5'); set(myfig.text.uL, 'string','0.5'); set(myfig.text.uR, 'string','-0.5'); set(myfig.text.xmax, 'string','10'); set(myfig.text.xmin, 'string','-10'); set(myfig.text.tfinal, 'string','5'); set(myfig.text.xprec, 'string','50'); set(myfig.text.Uic, 'string','blob(-5,-4)'); set(myfig.text.Sic, 'string','0'); set(myfig.text.topx,'string','x'); set(myfig.text.topy,'string','u_'); set(myfig.text.bottomx,'string','x'); set(myfig.text.bottomy,'string','s_'); set(myfig.menu,'Value',1); set(myfig.menu2,'Value',1); set(myfig.radio.top.auto, 'Value',1); set(myfig.radio.top.fixed, 'Value',0); set(myfig.radio.bottom.auto, 'Value',1); set(myfig.radio.bottom.fixed, 'Value',0); return; case 2 set(myfig.check.wave,'Value',1); set(myfig.check.pert,'Value',1); set(myfig.check.save,'Value',0); set(myfig.check.viscoterm,'Value',0); set(myfig.label.Uic,'string','e(x,t)'); set(myfig.label.Sic,'string','d(x,t)'); set(myfig.text.alpha, 'string','1') set(myfig.text.beta, 'string','5'); set(myfig.text.uL, 'string','0.5'); set(myfig.text.uR, 'string','-0.5'); set(myfig.text.xmax, 'string','10'); set(myfig.text.xmin, 'string','-10'); set(myfig.text.tfinal, 'string','5'); set(myfig.text.xprec, 'string','50'); set(myfig.text.Uic, 'string','blob(-6,-5)+blob(-5,-4)'); set(myfig.text.Sic, 'string','0'); set(myfig.text.topx,'string','x'); set(myfig.text.topy,'string','u_'); set(myfig.text.bottomx,'string','x'); set(myfig.text.bottomy,'string','s_'); set(myfig.menu,'Value',1); set(myfig.menu2,'Value',1); set(myfig.radio.top.auto, 'Value',1); set(myfig.radio.top.fixed, 'Value',0); set(myfig.radio.bottom.auto, 'Value',1); set(myfig.radio.bottom.fixed, 'Value',0); return; case 3 set(myfig.check.wave,'Value',1); set(myfig.check.pert,'Value',1); set(myfig.check.save,'Value',0); set(myfig.check.viscoterm,'Value',0); set(myfig.label.Uic,'string','e(x,t)'); set(myfig.label.Sic,'string','d(x,t)'); set(myfig.text.alpha, 'string','1') set(myfig.text.beta, 'string','5'); set(myfig.text.uL, 'string','0.5'); set(myfig.text.uR, 'string','-0.5'); set(myfig.text.xmax, 'string','10'); set(myfig.text.xmin, 'string','-10'); set(myfig.text.tfinal, 'string','3'); set(myfig.text.xprec, 'string','50'); set(myfig.text.Uic, 'string','0.5*sin(2*pi*x*(1))$-6,6$'); set(myfig.text.Sic, 'string','0'); set(myfig.text.topx,'string','x'); set(myfig.text.topy,'string','u_'); set(myfig.text.bottomx,'string','x'); set(myfig.text.bottomy,'string','s_'); set(myfig.menu,'Value',1); set(myfig.menu2,'Value',1); set(myfig.radio.top.auto, 'Value',1); set(myfig.radio.top.fixed, 'Value',0); set(myfig.radio.bottom.auto, 'Value',1); set(myfig.radio.bottom.fixed, 'Value',0); return; case 4 set(myfig.check.wave,'Value',1); set(myfig.check.pert,'Value',0); set(myfig.check.save,'Value',0); set(myfig.check.viscoterm,'Value',0); set(myfig.label.Uic,'string','e(x,t)'); set(myfig.label.Sic,'string','d(x,t)'); set(myfig.text.alpha, 'string','0.3') set(myfig.text.beta, 'string','5'); set(myfig.text.uL, 'string','0.5'); set(myfig.text.uR, 'string','-0.5'); set(myfig.text.xmax, 'string','2'); set(myfig.text.xmin, 'string','-2'); set(myfig.text.tfinal, 'string','4'); set(myfig.text.xprec, 'string','200'); set(myfig.text.Uic, 'string','waveu(0.5,-0.5,1,10,0)'); set(myfig.text.Sic, 'string','waves(0.5,-0.5,1,10,0)'); set(myfig.text.topx,'string','x'); set(myfig.text.topy,'string','u_'); set(myfig.text.bottomx,'string','x'); set(myfig.text.bottomy,'string','s_'); set(myfig.menu,'Value',2); set(myfig.menu2,'Value',1); set(myfig.radio.top.auto, 'Value',1); set(myfig.radio.top.fixed, 'Value',0); set(myfig.radio.bottom.auto, 'Value',1); set(myfig.radio.bottom.fixed, 'Value',0); return; case 5 set(myfig.check.wave,'Value',0); set(myfig.check.pert,'Value',1); set(myfig.check.save,'Value',0); set(myfig.check.viscoterm,'Value',0); set(myfig.label.Uic,'string','e(x,t)'); set(myfig.label.Sic,'string','d(x,t)'); set(myfig.text.alpha, 'string','1') set(myfig.text.beta, 'string','20'); set(myfig.text.uL, 'string','0.5'); set(myfig.text.uR, 'string','-0.5'); set(myfig.text.xmax, 'string','2'); set(myfig.text.xmin, 'string','-2'); set(myfig.text.tfinal, 'string','1'); set(myfig.text.xprec, 'string','200'); set(myfig.text.Uic, 'string','0'); set(myfig.text.Sic, 'string','0'); set(myfig.text.topx,'string','x'); set(myfig.text.topy,'string','u_'); set(myfig.text.bottomx,'string','x'); set(myfig.text.bottomy,'string','s_'); set(myfig.menu,'Value',1); set(myfig.menu2,'Value',1); set(myfig.radio.top.auto, 'Value',1); set(myfig.radio.top.fixed, 'Value',0); set(myfig.radio.bottom.auto, 'Value',1); set(myfig.radio.bottom.fixed, 'Value',0); return; case 6 set(myfig.check.wave,'Value',0); set(myfig.check.pert,'Value',1); set(myfig.check.save,'Value',0); set(myfig.check.viscoterm,'Value',0); set(myfig.label.Uic,'string','e(x,t)'); set(myfig.label.Sic,'string','d(x,t)'); set(myfig.text.alpha, 'string','0.4') set(myfig.text.beta, 'string','1'); set(myfig.text.uL, 'string','0.5'); set(myfig.text.uR, 'string','-0.5'); set(myfig.text.xmax, 'string','0.2'); set(myfig.text.xmin, 'string','-0.2'); set(myfig.text.tfinal, 'string','1'); set(myfig.text.xprec, 'string','1000'); set(myfig.text.Uic, 'string','0'); set(myfig.text.Sic, 'string','0'); set(myfig.text.topx,'string','x'); set(myfig.text.topy,'string','u_'); set(myfig.text.bottomx,'string','x'); set(myfig.text.bottomy,'string','s_'); set(myfig.menu,'Value',3); set(myfig.menu2,'Value',1); set(myfig.radio.top.auto, 'Value',1); set(myfig.radio.top.fixed, 'Value',0); set(myfig.radio.bottom.auto, 'Value',1); set(myfig.radio.bottom.fixed, 'Value',0); return; case 7 set(myfig.check.wave,'Value',0); set(myfig.check.pert,'Value',1); set(myfig.check.save,'Value',0); set(myfig.check.viscoterm,'Value',0); set(myfig.label.Uic,'string','e(x,t)'); set(myfig.label.Sic,'string','d(x,t)'); set(myfig.text.alpha, 'string','0.3') set(myfig.text.beta, 'string','1'); set(myfig.text.uL, 'string','0.5'); set(myfig.text.uR, 'string','-0.5'); set(myfig.text.xmax, 'string','0.1'); set(myfig.text.xmin, 'string','-0.1'); set(myfig.text.tfinal, 'string','1'); set(myfig.text.xprec, 'string','1000'); set(myfig.text.Uic, 'string','0'); set(myfig.text.Sic, 'string','0'); set(myfig.text.topx,'string','x'); set(myfig.text.topy,'string','u_'); set(myfig.text.bottomx,'string','x'); set(myfig.text.bottomy,'string','s_'); set(myfig.menu,'Value',1); set(myfig.menu2,'Value',1); set(myfig.radio.top.auto, 'Value',1); set(myfig.radio.top.fixed, 'Value',0); set(myfig.radio.bottom.auto, 'Value',1); set(myfig.radio.bottom.fixed, 'Value',0); return; end end % The functions below are used in the taylor approximation of the solution. % They take in the solution to both u and s at the previous time step, % along with the set of spatial indices over which the solution is defined % (i.e. everwhere except the ghost cells). Using the corresponding % relationships between the t-derivaties and x-derivatives, these functions % use centered finite difference approximations on the x-derivatives to % approximate the corresponding t-derrivatives. The dx and dt terms are % dealt with outside of these functions. % Relation: u_t = s_x - u*u_x function du = ut(s1,u12) du = s1 - u12; end % Relation: u_tt = (a^2 + s + u^2)*u_xx - 2*u*s_xx - u_x*s_x + 2*u*(u_x)^2 function du = utt(u0,u1,u2,s0,s1,s2,a) du = (a^2+s0+u0.^2).*u2 - 2*u0.*s2 - u1.*s1 + 2*u0.*u1.^2; end % Relation: s_t = (a^2 + s)*u_x - u*s_x function ds = st(u0,u1,s0,s1,a) ds = (a^2 + s0).*u1 - u0.*s1; end % Relation: s_tt = (a^2 + s + u^2)*s_xx - 2*u*(a^2 + s)*u_xx - (s_x)^2 function ds = stt(u0,u2,s0,s1,s2,a) ds = (a^2+s0+u0.^2).*s2 - 2*u0.*(a^2+s0).*u2 - s1.^2; end function du = tut(u12, v1, ebt) du = v1*ebt-u12; end function dv = vt(u0,v0,u1,v1,a,ebt) dv = (a^2*ebt+v0).*u1-u0.*v1; end function du = tutt(u0,v0,u1,v1,u2,v2,a,b,ebt) du = (a^2+v0/ebt+u0.^2).*u2-2*u0.*v2/ebt+2*u0.*u1.^2-v1.*u1/ebt-b*v1/ebt; end function dv = vtt(u0,v0,u1,v1,u2,v2,u12,u22,a,b,ebt) dv = -(a^2*ebt+v0).*u0.*u2 + (a^2+v0/ebt+u0.^2).*v2 - u0.*v1.*u1 + ... (a^2*ebt + v0).*u1.^2 - v1.^2/ebt - (a^2*ebt+v0).*u22 + v1.*u12 + ... b*a^2*ebt*u1; end % The function below checks to make sure that all of the parameters the % user has entered are valid. If something is incorrect the a warning % message is displayed and the program is haluted. The program runs a basic % check on the user specified initial conditions, but they are checked more % rigorously later on. function valid = validconditions(fig) msg1 = 'Bad parameter settings!'; msg2 = 'Please specify initial conditions'; msg3 = 'xmin must be less than xmax!'; msg4 = 'tend must be greater than zero!'; msg5 = 'beta must be greater than zero!'; msg6 = '$,$ mismatch in initial conditions'; xmin = get(fig.text.xmin,'string'); xmax = get(fig.text.xmax,'string'); tend = get(fig.text.tfinal,'string'); alpha = get(fig.text.alpha,'string'); beta = get(fig.text.beta,'string'); uL = get(fig.text.uL,'string'); uR = get(fig.text.uR,'string'); U = get(fig.text.Uic,'string'); S = get(fig.text.Sic,'string'); if (isempty(xmin) | isempty(xmax) | isempty(tend) | isempty(alpha)|isempty(beta) | ... isempty(uL) | isempty(uR)) errmsg(msg1); valid = false; return; end if (isempty(U) | isempty(S)) errmsg(msg2); valid = false; return; end if (str2num(xmin) >= str2num(xmax)) errmsg(msg3); valid = false; return; end if (str2num(tend) <= 0) errmsg(msg4); valid = false; return; end if (beta <= 0) errmsg(msg5); valid = false; return; end if (isempty(parse(U,S))) errmsg(msg6); valid = false; return; end valid = true; end % This function simply displays a warning message in the status box where % the message is specified by the input string. function errmsg(string) global myfig; set(myfig.status,'string',string,'ForegroundColor',[1 0 0]); getframe(gcf); end function status(string) global myfig; set(myfig.status,'string',string,'ForegroundColor',[0 0 0]); getframe(gcf); end function progress(string) global myfig; set(myfig.progress,'string',string,'ForegroundColor',[0 0 0]); getframe(gcf); end % The parse function takes as input the two strings specifying the u and s % initial conditions entered by the user. It clears up all white spaces % and makes sure all operations are vector operations. The main function % is to input the piecewise solution into a structure so that it can be % more easily calculated in the parent function. The user may enter a % piecewise function as follows f(x) = expr1.. $-1,1$ expr2.. $1,5$ where % the valid interval of the particular piece of the piecewise function is % specified within the dollar signs. If the limits were not entered % correctly the function displays and error message and returns an empty % set of initial condtions. function ic = parse(uic,sic) uic = stringPrep(uic); sic = stringPrep(sic); [trash n1] = size(findstr(uic,'$')); [trash n2] = size(findstr(sic,'$')); m1 = n1/2; m2 = n2/2; if (m1 ~= floor(m1) | m2 ~= floor(m2)) errmsg('Invalid initial conditions'); ic = struct([]); return; end ic.u.size = m1; ic.s.size = m2; if (m1 ~= 0) rem = uic; for i = 1:m1 [ic.u.exp(i).str rem] = strtok(rem,'$'); rem = rem(2:length(rem)); [llim rem] = strtok(rem, ','); rem = rem(2:length(rem)); [rlim rem] = strtok(rem,'$'); rem = rem(2:length(rem)); ic.u.llim(i) = str2num(llim); ic.u.rlim(i) = str2num(rlim); if (isempty(ic.u.llim(i)) | isempty(ic.u.rlim(i))) errmsg('Invalid initial conditions!'); ic = struct([]); return; end end else ic.u.exp(1).str = uic; end if (m2 ~= 0) rem = sic; for i = 1:m2 [ic.s.exp(i).str rem] = strtok(rem,'$'); rem = rem(2:length(rem)); [llim rem] = strtok(rem, ','); rem = rem(2:length(rem)); [rlim rem] = strtok(rem,'$'); rem = rem(2:length(rem)); ic.s.llim(i) = str2num(llim); ic.s.rlim(i) = str2num(rlim); if (isempty(ic.s.llim(i)) | isempty(ic.s.rlim(i))) errmsg('Invalid initial conditions!'); ic = struct([]); return; end end else ic.s.exp(1).str = sic; end end function y = blob(a,b) global grid; x = grid.x; if (abs(a-b) < 1 & a ~= b) y = (x >= a & x <= b).*sin((x-a)*pi/(a-b)); end if (abs(a-b) >=1) y = (x >=a & x < a+0.5).*(0.5-0.5*cos(2*pi*(x-a))) + ... (x >= a+0.5 & x < b-0.5).*(1) + ... (x >= b-0.5 & x < b).*(0.5+0.5*cos(2*pi*(x-(b-0.5)))); end if (a == b) y = zeros(size(x)); end end function u = waveu(uL,uR,a,b,center) global grid; global parameters; if (center < 0) xInterval = [grid.xmin grid.xmax-center]; else xInterval = [grid.xmin-center grid.xmax]; end [s u] = wavesolution(a,b,uL,uR,xInterval,grid.x-center); end function s = waves(uL,uR,a,b,center) global grid; global parameters; if (center < 0) xInterval = [grid.xmin grid.xmax-center]; else xInterval = [grid.xmin-center grid.xmax]; end [s u] = wavesolution(a,b,uL,uR,xInterval,grid.x-center); end function str2 = stringPrep(str1) str2 = strrep(str1,' ',''); str2 = strrep(str2,'*','.*'); str2 = strrep(str2,'/','./'); str2 = strrep(str2,'^','.^'); str2 = strrep(str2,'..','.'); end function str2 = stringPrepAx(str1,bool) str2 = stringPrep(str1); switch (bool) case false str2 = strrep(str2,'s_','sol.s'); str2 = strrep(str2,'u_','sol.u'); str2 = strrep(str2,'sw_','sol.sw'); str2 = strrep(str2,'uw_','sol.uw'); str2 = strrep(str2,'e_','(sol.u-sol.uw)'); str2 = strrep(str2,'d_','(sol.s - sol.sw)'); str2 = strrep(str2,'Lw','(sol.uw-2*sqrt(sol.sw+a^2))'); str2 = strrep(str2,'Rw','(sol.uw+2*sqrt(sol.sw+a^2))'); str2 = strrep(str2,'L','(sol.u-2*sqrt(sol.s+a^2))'); str2 = strrep(str2,'R','(sol.u+2*sqrt(sol.s+a^2))'); return; case true str2 = strrep(str2,'s_','sol.sw'); str2 = strrep(str2,'u_','sol.uw'); str2 = strrep(str2,'e_','sol.z'); str2 = strrep(str2,'d_','sol.z'); str2 = strrep(str2,'Lw','(sol.uw-2*sqrt(sol.sw+a^2))'); str2 = strrep(str2,'Rw','(sol.uw+2*sqrt(sol.sw+a^2))'); str2 = strrep(str2,'L','(sol.uw-2*sqrt(sol.sw+a^2))'); str2 = strrep(str2,'R','(sol.uw+2*sqrt(sol.sw+a^2))'); end end function str2 = stringPrepAx2(str1,bool) str2 = stringPrep(str1); switch (bool) case false str2 = strrep(str2,'s_','sol2.s'); str2 = strrep(str2,'u_','sol2.u'); str2 = strrep(str2,'e_','(sol2.u-sol2.uw)'); str2 = strrep(str2,'d_','(sol2.s - sol2.sw)'); str2 = strrep(str2,'Lw','(sol2.uw-2*sqrt(sol2.sw+a^2))'); str2 = strrep(str2,'Rw','(sol2.uw+2*sqrt(sol2.sw+a^2))'); str2 = strrep(str2,'L','(sol2.u-2*sqrt(sol2.s+a^2))'); str2 = strrep(str2,'R','(sol2.u+2*sqrt(sol2.s+a^2))'); return; case true str2 = strrep(str2,'s_','sol2.sw'); str2 = strrep(str2,'u_','sol2.uw'); str2 = strrep(str2,'e_','sol2.z'); str2 = strrep(str2,'d_','sol2.z'); str2 = strrep(str2,'Lw','(sol2.uw-2*sqrt(sol2.sw+a^2))'); str2 = strrep(str2,'Rw','(sol2.uw+2*sqrt(sol2.sw+a^2))'); str2 = strrep(str2,'L','(sol2.uw-2*sqrt(sol2.sw+a^2))'); str2 = strrep(str2,'R','(sol2.uw+2*sqrt(sol2.sw+a^2))'); end end function y = dx(m) global grid; global playopts; [a b] = size(m); i = 2:a-1; n = 1:b; y(i,n) = (m(i+1,n) - m(i-1,n))/(2*grid.dx); y(1,n) = (m(1+1,n) - m(1,n))/(grid.dx); y(a,n) = (m(a,n) - m(a-1,n))/(grid.dx); end function y = dt(m) global grid; global playopts; [a b] = size(m); i = 1:a; n = 2:b-1; y(i,n) = (m(i,n+1) - m(i,n-1))/(index2time(n+1) - index2time(n-1)); y(i,1) = (m(i,1+1) - m(i,1))/(index2time(1+1)-index2time(1)); y(i,b) = (m(i,b) - m(i,b-1))/(index2time(b)-index2time(b-1)); end function y = dxx(m) global gird; global playopts; [a b] = size(m); i = 2:a-1; n = 1:b; y(i,n) = (m(i+1,n) - 2*m(i,n) + m(i-1,n))/(grid.dx^2); y(1,n) = (2*m(1,n) - 5*m(1+1,n) + 4*m(1+2,n) - m(1+3,n))/(grid.dx^2); y(a,n) = (-m(a-3,n) + 4*m(a-2,n) - 5*m(a-1,n) + 2*m(a,n))/(grid.dx^2); end function index = time2index(time) global grid; global playopts; for i = 1:grid.T if (abs(index2time(i) - time) <= 1/playopts.fps) index = i; return; end end index = grid.T; end function time = index2time(index) global grid; global playopts; k = 1 + floor(grid.tprec/playopts.fps*(index-1)); time = k*grid.dt; end function index = x2index(x) global grid; index = floor(grid.xprec*(x - grid.xmin))+1; if (index < grid.x0) index = 1; end if (index > grid.xN) index = grid.XN; end end function x = index2x(index) global grid; x = grid.xmin + (index-1)/grid.xprec; end function y = int(m) global grid; y = sum(m)*grid.dx; end function y = intt(m,a,b) global grid; indexA = time2index(a); indexB = time2index(b); t = indexA:indexB; y = sum(m(:,t))*grid.dt; end function y = energy(u,s) global grid; global parameters; y = u.^2/2 - (parameters.uL + parameters.uR)/2*u - s; end