%%%%%%%%%%%%%%%%%%%% Source code from this run of demo1.m: %%%%%%%%%%%%%%%%%%%% %% NOTE: Nt = 5600 timesteps! %% % Physical parameters domain_width = 1; % domain size c = 1; % wave speed tmax = 35; % total time for simulation % Grid setup Nx = 32; Delta = domain_width/Nx; x = Delta*(0:(Nx-1)) - domain_width/2; % Timestep setup SAFETY = 5 % factor by which to exceed the courant condition, h < Delta/c. hmax = Delta/(c*SAFETY); % biggest timestep that would satisfy SAFETY. Nt = tmax * ceil(1/hmax); % number of timesteps to reach tmax, with h < hmax, % ensuring that all integer values of t are on the time grid. h = tmax/Nt; % resulting timestep, h times = linspace(h, tmax, Nt); % timestep array, leaving out the initial t=0. h_dimensionless = h * c / Delta % dimensionless timestep epsilon = h_dimensionless^2 + 4*h_dimensionless^4 % parameter required by rk2plus. Rule of thumb based on an approximation % of the stability diagram. % Initial conditions pulse_width = domain_width/16; % std dev of Gaussian pulse. u0 = exp(-x.^2/(2*pulse_width^2)); % initial wavefunction. u0_x = - x.*u0/(pulse_width^2); % partial derivative of u0 wrt x. udot0 = - u0_x * c; % pulse propagates to the right. q0 = [u0, udot0]; % entire initial state vector. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Output: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% octave:1> demo1 SAFETY = 5 h_dimensionless = 0.20000 epsilon = 0.046400 cputime_elapsed_min_rk2plus = 11.647 cputime_elapsed_min_rk2 = 12.410