ssa_coeffs.mac (2535B)
1 /* -*- mode: maxima -*- */ 2 isolate_wrt_times:true$ 3 4 /* Shifts in x and y directions. */ 5 shift_x(expr, d) := op(expr)[args(expr)[1]+d, args(expr)[2]]$ 6 shift_y(expr, d) := op(expr)[args(expr)[1], args(expr)[2]+d]$ 7 shift[x] : shift_x; shift[y] : shift_y; 8 9 /* weights */ 10 weight(var, direction, d) := shift[direction](subst(w, op(var), var), d)$ 11 12 /* one-sided differences with weights */ 13 d_px(var) := weight(var, x, 1/2) * (shift[x](var, 1) - var)$ 14 d_mx(var) := weight(var, x, -1/2) * (var - shift[x](var, -1))$ 15 16 d_py(var) := weight(var, y, 1/2) * (shift[y](var, 1) - var)$ 17 d_my(var) := weight(var, y, -1/2) * (var - shift[y](var, -1))$ 18 19 /* centered differences defined as sums of one-sided differences with weights */ 20 D_x(foo) := d_px(foo) + d_mx(foo)$ 21 D_y(foo) := d_py(foo) + d_my(foo)$ 22 23 load("ssa.mac")$ 24 25 denominator : (4 * dx**2 * dy**2)$ 26 27 /* Clear the denominator */ 28 e1 : ''lhs1 * denominator, expand; 29 e2 : ''lhs2 * denominator, expand; 30 31 /* define weights to be equal to 1 in the interior; give them PIK names otherwise */ 32 a_ones : [w[i+1/2,j] = 1, 33 w[i-1/2,j] = 1, 34 w[i,j+1/2] = 1, 35 w[i,j-1/2] = 1, 36 N[i+1/2,j] = c_e, 37 N[i-1/2,j] = c_w, 38 N[i,j+1/2] = c_n, 39 N[i,j-1/2] = c_s, 40 w[i-1/2,j+1] = 1, 41 w[i+1/2,j+1] = 1, 42 w[i+1,j+1/2] = 1, 43 w[i+1,j-1/2] = 1, 44 w[i+1/2,j-1] = 1, 45 w[i-1/2,j-1] = 1, 46 w[i-1,j-1/2] = 1, 47 w[i-1,j+1/2] = 1]$ 48 49 a_names : [w[i+1/2,j] = aPP, 50 w[i-1/2,j] = aMM, 51 w[i,j+1/2] = bPP, 52 w[i,j-1/2] = bMM, 53 N[i+1/2,j] = c_e, 54 N[i-1/2,j] = c_w, 55 N[i,j+1/2] = c_n, 56 N[i,j-1/2] = c_s, 57 w[i-1/2,j+1] = aMn, 58 w[i+1/2,j+1] = aPn, 59 w[i+1,j+1/2] = bPe, 60 w[i+1,j-1/2] = bMe, 61 w[i+1/2,j-1] = aPs, 62 w[i-1/2,j-1] = aMs, 63 w[i-1,j-1/2] = bMw, 64 w[i-1,j+1/2] = bPw]$ 65 66 if interior then (e1 : at(e1, a_ones), e2 : at(e2, a_ones)) 67 else (e1 : at(e1, a_names), e2 : at(e2, a_names))$ 68 69 /* Finally compute coefficients: */ 70 for m: -1 thru 1 do (for n: -1 thru 1 do (c1u[m,n] : combine(expand(coeff(e1, u[i+m,j+n]) / denominator)))); 71 for m: -1 thru 1 do (for n: -1 thru 1 do (c1v[m,n] : combine(expand(coeff(e1, v[i+m,j+n]) / denominator)))); 72 73 for m: -1 thru 1 do (for n: -1 thru 1 do (c2u[m,n] : combine(expand(coeff(e2, u[i+m,j+n]) / denominator)))); 74 for m: -1 thru 1 do (for n: -1 thru 1 do (c2v[m,n] : combine(expand(coeff(e2, v[i+m,j+n]) / denominator)))); 75