clear clc clg echo on % ------------------------- WaveEqn.m ------------------------- % % The Wave Equation and Its Numerical Solution % Lecture 4 % % % % % % Chem 455 % Summer Quarter 1997 % Instructor - J. B. Callis % % Press any key to continue pause clc % First establish the constants, initial conditions % and execution parameters: % % For the initial condition, imagine a string L = 100 cm long, % that is stretched 0.1 cm in the middle. L = 100; Dx = 5; % Spatial step in cm x = 0:Dx:L; % This is the vector of spatial coordinate values % % % Press any key to continue pause clc % Now use the if-else-end structure to calculate the initial % values for i = 1:length(x) % The for loop runs over all values of x if x(i)<=L/2 % Amplitude is increasing over this interval y(i)=0.1*x(i)/(L/2); else y(i) = 0.1*(L-x(i))/(L/2); % Amplitude is decreasing end end % % % Press any key to see a plot of the initial condition pause clc plot(x,y) title('Initial Condition for Vibrating String') xlabel('Distance Along String, cm') ylabel('Amount String is Initially Displaced, cm') pause clc % Suppose the string is stretched until the wave velocity is % 5e4 cm/sec. This fixes the relationship between the spatial % and temporal steps according to: % % c = Dx/Dt, or In this case since the spatial steps are 5 cm, % time steps are Dt = Dx/c % c = 5e4; Dt =Dx/c; % Calculate the time steps % % % Press any key to continue pause clc % How far shall we take the time? Certainly, we want to go % over one cycle of the lowest frequency of the vibration. % This period may be calculated from time = 2*L/c % sim_time = 2*L/c; Nsteps = round(sim_time/Dt); % round to nearest integer u = zeros(Nsteps,length(x)); % Set up (pre-allocate) the u matrix % % Press any key to continue pause clc % We now are ready to propagate our difference formula from % the initial condition, u(x,0): % % the first time step is special since we do not know the value % of u(x,-Dt). However, since the initial condition represents % the maximum extent of displacement, and since, the displacement % oscillates, we can use the equality: u(x,-Dt) = u(x,Dt) on the % first time step only. % % In addition we must recall the boundary conditions that % the string displacement must be 0 outside the interval % 0<=x<=L. % % This will require some logical statements. % Press any key to continue pause clc % First we calculate the initial string displacement in row 1: u(1,:) = y; % Put y in the first row of the u matrix % now calculate the second row for i = 2:(length(x)-1) % boundary positions, x=0, L already zero u(2,i)=(1/2)*(u(1,i+1)+u(1,i-1)); % calculation for first % step only end % calculate the rest of the rows for j = 2:Nsteps for i = 2:(length(x)-1) u(j+1,i)=u(j,i+1)+u(j,i-1)-u(j-1,i); end end % Calculation is finished. Press any key to play a movie pause clg for j = 1:Nsteps plot(x,u(j,:)) axis([0 max(x) -.1 .1]); pause(.1) end pause %Press any key to see a plot of all of the data clg plot(x,u) title('Amplitude of String as a Function of Time') xlabel('Distance Along String, cm') ylabel('Amount String is Displaced, cm') pause