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 multi-step methods The third-order Adams Moulton method uses the time-levels t n +1 , t n , and t n − 1 to interpolate the right hand side function f to third-order accuracy using quadratic polynomials. Derive the coefficients a i of this third-order 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 predictor-corrector approach where the predicted value is obtained from a third-order explicit scheme and used in the evaluation of the right hand side term f n +1 . 2 Leap-Frog Trapezoidal time stepping One remedy to the computational mode of the Leap-Frog 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 leap-frog 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 Leap-Frog, 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=(xmax-xmin)/(M-1); x =xmin:dx:xmax; N=201; ymin = 0; ymax=2.; dy=(ymax-ymin)/(N-1); 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.^2-4*a.*c); A1 = (-b+det)./(2*a); A2 = (-b-det)./(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 non-dimensional 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
- Trigraph, Subroutine, right hand, r8, Right-hand rule