Attributes[buildmatrix] = {HoldRest} buildmatrix::usage = "buildmatrix[\"f\",m] builds a matrix m from the file f \n with the structure:\n dimension \n row index column index value \n ..." buildmatrix/: buildmatrix[file_, matrix_] := Block[{data, f, i, j, l, n}, f = OpenRead[file]; n = Read[f, Number]; matrix = DiagonalMatrix[Table[0,{n}]]; data = ReadList[f, {Number, Number, Number}]; Do[i = data[[l,1]]; j = data[[l,2]]; matrix[[i,j]] = data[[l,3]], {l, Length[data]}]; Close[f]] Attributes[buildpat] = {HoldRest} buildpat::usage = "buildpat[\"f\",m] builds a matrix pattern m from the file f \n with the structure:\n dimension \n row index column index value \n ..." buildpat/: buildpat[file_, matrix_] := Block[{data, f, i, j, l, n}, f = OpenRead[file]; n = Read[f, Number]; matrix = DiagonalMatrix[Table[0,{n}]]; data = ReadList[f, {Number, Number, Number}]; Do[i = data[[l,1]]; j = data[[l,2]]; matrix[[i,j]] = 1.0, {l, Length[data]}]; Close[f]; matrix = Reverse[ Map[ If[# != 0, 0, 1]&, matrix, {2}] ]; ] lcp/: lcp[m_, q_, n_] := Block[{l, indices, t, x}, indices = Flatten[Apply[Outer, Prepend[Table[{l, l + n}, {l, n}], List]], n-1]; t = Table[Join[IdentityMatrix[n][[l]],-m[[l]] ], {l,n}]; Do[x = Inverse[Transpose[Transpose[t][[indices[[l]] ]] ] ] . q; Print[indices[[l]], " ", x]; Null, {l, 2^n}] ] lcp::usage = "lcp prints all the basic sequences and basic variables for the n x n nondegenerate LCP (m/q)." (* Formulas on 21 point stencil. *) v[x_,y_] = Normal[Series[u[x,y],{x,0,5},{y,0,5}]]; sumx= Expand[c1 v[-1,2] +c2 v[0,2] +c3 v[1,2] + c4 v[-2,1] + c5 v[-1,1] + c6 v[0,1] + c7 v[1,1] + c8 v[2,1] + c9 v[-2,0] + c10 v[-1,0] + c11 v[0,0] + c12 v[1,0] + c13 v[2,0] + c14 v[-2,-1] + c15 v[-1,-1] + c16 v[0,-1] + c17 v[1,-1] + c18 v[2,-1] + c19 v[-1,-2] + c20 v[0,-2] + c21 v[1,-2] ] ; eq[ 1]=Coefficient[sumx,Derivative[0,0][u][0,0]] == 0; eq[ 2]=Coefficient[sumx,Derivative[1,0][u][0,0]] == 0; eq[ 3]=Coefficient[sumx,Derivative[0,1][u][0,0]] == 0; eq[ 4]=Coefficient[sumx,Derivative[2,0][u][0,0]] == 0; eq[ 5]=Coefficient[sumx,Derivative[1,1][u][0,0]] == 0; eq[ 6]=Coefficient[sumx,Derivative[0,2][u][0,0]] == 0; eq[ 7]=Coefficient[sumx,Derivative[3,0][u][0,0]] == 0; eq[ 8]=Coefficient[sumx,Derivative[2,1][u][0,0]] == 0; eq[ 9]=Coefficient[sumx,Derivative[1,2][u][0,0]] == 0; eq[10]=Coefficient[sumx,Derivative[0,3][u][0,0]] == 0; eq[11]=Coefficient[sumx,Derivative[4,0][u][0,0]] == 0; eq[12]=Coefficient[sumx,Derivative[3,1][u][0,0]] == 0; eq[13]=Coefficient[sumx,Derivative[2,2][u][0,0]] == 0; eq[14]=Coefficient[sumx,Derivative[1,3][u][0,0]] == 1; eq[15]=Coefficient[sumx,Derivative[0,4][u][0,0]] == 0; eq[16]=Coefficient[sumx,Derivative[5,0][u][0,0]] == 0; eq[17]=Coefficient[sumx,Derivative[4,1][u][0,0]] == 0; eq[18]=Coefficient[sumx,Derivative[3,2][u][0,0]] == 0; eq[19]=Coefficient[sumx,Derivative[2,3][u][0,0]] == 0; eq[20]=Coefficient[sumx,Derivative[1,4][u][0,0]] == 0; eq[21]=Coefficient[sumx,Derivative[0,5][u][0,0]] == 0; Solve[Prepend[Table[eq[i],{i,21}],c3==0]] (*Set some c's to eliminate some of the free parameters.*) (* Differential operator formulas for equation transformed from equilateral (x-y coordinates) triangle to right triangle (s-t coordinates). *) s/: s[x_, y_] := x - (1*y)/Sqrt[3] t/: t[x_, y_] := (2*y)/Sqrt[3] u/: u[x_, y_] := v[s[x, y], t[x, y]] ux/: ux[x_, y_] = Derivative[1, 0][v][x - y/3^(1/2), (2*y)/3^(1/2)] uy/: uy[x_, y_] = (2*Derivative[0, 1][v][x - y/3^(1/2), (2*y)/3^(1/2)])/3^(1/2) - Derivative[1, 0][v][x - y/3^(1/2), (2*y)/3^(1/2)]/3^(1/2) xytost = {x -> s + t/2, y -> (3^(1/2)*t)/2, v[s, t] -> vV[0, 0], Derivative[0, 1][v][s, t] -> vV[0, 1], Derivative[0, 2][v][s, t] -> vV[0, 2], Derivative[0, 3][v][s, t] -> vV[0, 3], Derivative[0, 4][v][s, t] -> vV[0, 4], Derivative[1, 0][v][s, t] -> vV[1, 0], Derivative[1, 1][v][s, t] -> vV[1, 1], Derivative[1, 2][v][s, t] -> vV[1, 2], Derivative[1, 3][v][s, t] -> vV[1, 3], Derivative[1, 4][v][s, t] -> vV[1, 4], Derivative[2, 0][v][s, t] -> vV[2, 0], Derivative[2, 1][v][s, t] -> vV[2, 1], Derivative[2, 2][v][s, t] -> vV[2, 2], Derivative[2, 3][v][s, t] -> vV[2, 3], Derivative[2, 4][v][s, t] -> vV[2, 4], Derivative[3, 0][v][s, t] -> vV[3, 0], Derivative[3, 1][v][s, t] -> vV[3, 1], Derivative[3, 2][v][s, t] -> vV[3, 2], Derivative[3, 3][v][s, t] -> vV[3, 3], Derivative[3, 4][v][s, t] -> vV[3, 4], Derivative[4, 0][v][s, t] -> vV[4, 0], Derivative[4, 1][v][s, t] -> vV[4, 1], Derivative[4, 2][v][s, t] -> vV[4, 2], Derivative[4, 3][v][s, t] -> vV[4, 3], Derivative[4, 4][v][s, t] -> vV[4, 4]} del2/: del2[f_] := Through[((D[#1, {x, 2}] & ) + (D[#1, {y, 2}] & ))[f], Plus] (* Laplacian operator with respect to x, y, applied to function f. *) (* Differential operator formulas for equation transformed from trapezoid (x-y coordinates) to unit square (s-t coordinates). *) s/: s[x_, y_] := (d y + 3 x - 3 d)/(2(d y + 3 Sqrt[3] - 3 d)) t/: t[x_, y_] := y/3 u/: u[x_, y_] := v[s[x, y], t[x, y]] ux/: ux[x_, y_] = (Derivative[1, 0][v][ (d y + 3 x - 3 d)/(2(d y + 3 Sqrt[3] - 3 d)), y/3] D[s[x,y], x]) uy/: uy[x_, y_] = (Derivative[1, 0][v][ (d y + 3 x - 3 d)/(2(d y + 3 Sqrt[3] - 3 d)), y/3] D[s[x,y], y]) + (1/3) (Derivative[0, 1][v][ (d y + 3 x - 3 d)/(2(d y + 3 Sqrt[3] - 3 d)), y/3] ) xytost = Flatten[{ (d y + 3 x - 3 d)/(2(d y + 3 Sqrt[3] - 3 d)) -> ss, x -> d + 2 Sqrt[3] ss - 2 d ss - d tt + 2 d ss tt, y -> 3 tt, Table[Derivative[i, j][v][ss, tt] -> vV[i, j], {i, 0, 4}, {j, 0, 4}] }] del2/: del2[f_] := Through[((D[#1, {x, 2}] & ) + (D[#1, {y, 2}] & ))[f], Plus] (* Laplacian operator with respect to x, y, applied to function f. *) (* Compute and save coefficients of partials in full differential operator. LU = del2[del2[u[x,y]]] + aA del2[ux[x,y]] + bB del2[uy[x,y]] + cC ux[x,y] + dD uy[x,y] //. xytost ; Do[Do[ cvV[i,j] = Simplify[Coefficient[Collect[ LU, vV[i,j]], vV[i,j]]]; OutputForm[SequenceForm["cvV[",i,",",j,"]=", FortranForm[cvV[i,j]]]] >>> coef.f , {j, 0, 4-i, 1}], {i, 0, 4, 1}]; *)