This preview has intentionally blurred sections. Sign up to view the full version.
View Full Document
Unformatted text preview: MPO 662 – Problem Set 3 1 Implicit multistep methods The thirdorder Adams Moulton method uses the timelevels t n +1 , t n , and t n − 1 to interpolate the right hand side function f to thirdorder accuracy using quadratic polynomials. Derive the coefficients a i of this thirdorder integration method: u n +1 = u n + Δ t parenleftBig a − 1 f n +1 + a f n + a 1 f n − 1 parenrightBig (1) Note that a modified version of this algorithm uses a predictorcorrector approach where the predicted value is obtained from a thirdorder explicit scheme and used in the evaluation of the right hand side term f n +1 . 2 LeapFrog Trapezoidal time stepping One remedy to the computational mode of the LeapFrog scheme is to combine it with the trapezoidal step. The combined scheme can be written as follows: u ∗ = u n − 1 + 2Δ tf ( u n , t n ) (2) u n +1 = u n + Δ tf parenleftbigg u n + u ∗ 2 , t n + 1 2 parenrightbigg (3) The scheme can be viewed as a predictor corrector method. The first step predicts u at time level n + 1 using a leapfrog scheme, while the second one is a trapezoidal step centered at n + 1 2 . Unlike the pure trapezoidal scheme, this one is explicit and uses the average of u n and u ∗ as an estimate of u n + 1 2 . Investigate the stability of the scheme for the case f = κu by deriving the expression for the amplification factor. Is there a parasitic mode? and if yes which one? Plot the amplification factor for κ Δ t = z = x + iy where  x  < 1 and 0 ≤ y ≤ 2. Compare its stability region to that of the LeapFrog, RK2 and AB2 schemes. Plotting stability diagrams is rather cumbersome especially if they involve the square roots of complex numbers. Here is a sample matlab script to achieve that. M=201; xmin =1; xmax=1.; dx=(xmaxxmin)/(M1); x =xmin:dx:xmax; N=201; ymin = 0; ymax=2.; dy=(ymaxymin)/(N1); y =ymin:dy:ymax; z= x’*ones(size(y)) + i * ones(size(x))’*y; % roots of equation a A^2 + b*A +c, where a,b and c are functions of z % a=; % b=; % c=; det = sqrt(b.^24*a.*c); A1 = (b+det)./(2*a); A2 = (bdet)./(2*a); A = max(abs(A1),abs(A2)); contour(x,y,A’,[1 1]); % will draw the contour of neutral stability 3 Linear shallow water equations The linearized nondimensional shallow water equations in a channel of unit depth and unit width are u t + η x = 0 (4) v t + η y = 0 (5) η t + u x + v y = 0 (6) ( u, v, η ) are the velocity, in the ( x, y ) direction and the hydrostatic pressure due to the motion of the free surface. the solution in space can be described by a Fourier series whose expansion takes the following form: u ( x, y, t ) = ˆ u ( t ) sin px cos qy (7) v ( x, y, t ) = ˆ v ( t ) cos px sin qy (8) η ( x, y, t ) = ˆ η ( t ) cos px cos qy (9) where ( p, q ) are wavenumbers in the (...
View
Full Document
 Spring '08
 Iskandarani,M
 Trigraph, Subroutine, right hand, r8, Righthand rule

Click to edit the document details