要素行列作成で利用する

弱形式:

行列表現:

Aとfの要素ごとの寄与分（要素行列）計算バッチファイル:

kill(all);
hk:sqrt(2);
B:genmatrix(lambda([i,j],b[i-1]),3,1);
C:genmatrix(lambda([i,j],c[i-1]),3,1);
G:(s/12)*matrix([2,1,1],[1,2,1],[1,1,2]);
 
M11:B.transpose(B);
M12:C.transpose(B);
M22:C.transpose(C);
S11:2*M11+M22;
S12:M12;
S22:M11+2*M22;
S:s*addrow(addcol(S11,S12),addcol(transpose(S12),S22));
 
H1:B.genmatrix(lambda([i,j],1),1,3);
H2:C.genmatrix(lambda([i,j],1),1,3);
H:-s/3*addrow(H1,H2);
 
J11:B.genmatrix(lambda([i,j],1),1,3);
J12:1/2*C.genmatrix(lambda([i,j],1),1,3);
J13:1/2*C.genmatrix(lambda([i,j],1),1,3);
J14:zeromatrix(3,3);
J21:zeromatrix(3,3);
J22:1/2*B.genmatrix(lambda([i,j],1),1,3);
J23:1/2*B.genmatrix(lambda([i,j],1),1,3);
J24:C.genmatrix(lambda([i,j],1),1,3);
J:s/3*addrow(addcol(addcol(addcol(J11,J12),J13),J14),addcol(addcol(addcol(J21,J22),J23),J24));
 
L2:addrow(addcol(G,zeromatrix(3,3)),addcol(zeromatrix(3,3),G));
L4:addrow(addcol(L2,zeromatrix(6,6)),addcol(zeromatrix(6,6),L2));
 
A11:etas*S+rho/tau*L2;
A12:H;
A13:J;
A21:transpose(H);
A22:alpha*hk^2/2/etap*s*(B.transpose(B)+C.transpose(C));
A23:zeromatrix(3,12);
A31:-2*etap*transpose(J);
A32:zeromatrix(12,3);
A33:(1+lambda/tau)*L4;
Ak:addrow(addrow(addcol(addcol(A11,A12),A13),addcol(addcol(A21,A22),A23)),addcol(addcol(A31,A32),A33));
 
uu:genmatrix(lambda([i,j],if i<=9 then rho/tau*uk[i-1] else lambda/tau*uk[i-1]),21,1);
ff:genmatrix(lambda([i,j],f[i-1]),21,1);
L6:addrow(addrow(addcol(L2,zeromatrix(6,15)),zeromatrix(3,21)),addcol(zeromatrix(12,9),L4));
fk:L6.(ff+uu);
 
load("cform.lisp");
with_stdout("Ak.txt",cform(Ak));
with_stdout("fk.txt",cform(fk));
 

バッチを処理した後Ak.txtとfk.txtの中身をみると行列のインデックスが1から始まってしまっています（フォートラン風に）
0からに変更するためにcform.lispを編集してもう一度バッチ処理します。ピンク色のI 0とJ 0が初期状態で1になっているので0に変更します

・
・
・
;; Takes a name and a matrix and prints a sequence of C assignment
;; statements of the form
;;  NAME[I][J] = <corresponding matrix element>
;;  The indcies I, J will be counted from 1.
 
(DEFMFUN $CFORMMX (NAME MAT &OPTIONAL (STREAM #-CL NIL #+CL *standard-output*)
			 &AUX ($LOADPRINT NIL) (K 'array))
  (COND ((NOT (symbolp NAME))
	 (MERROR "~%First argument to CFORMMX must be a symbol."))
	((NOT ($MATRIXP MAT))
	 (MERROR "Second argument to CFORMMX not a matrix: ~M" MAT)))
  (DO ((MAT (CDR MAT) (CDR MAT)) (I 0 (f1+ I))) ((NULL MAT))
      (DECLARE (FIXNUM I))
      (DO ((M (CDAR MAT) (CDR M)) (J 0 (f1+ J))) ((NULL M))
	  (DECLARE (FIXNUM J))
	  (C-PRINT `((MEQUAL) ((((,NAME ,K) ,I) ,K) ,J) ,(CAR M)) STREAM)))
  '$DONE)
 
;; End:
 

またfkはベクトルのつもりで作ったがmaximaでは内部的にベクトルも行列として処理するので、余計なインデックスがついてしまいます
そこでfk.txtをメモ帳等で開いて"][0]"を"]"に置換します
最後にAk.txtとfk.txtの内容をヘッダファイルにコピペします（非常に長くなるのでメインプログラムにそのままコピペするのはお勧めできません）

void make_Ak(double **Ak,double s,double b[],double c[],double tau){
Ak[0][0] = rho*s*pow(tau,-1)/6.0+(pow(c[0],2)+2*pow(b[0],2))*etas*s;
Ak[0][1] = rho*s*pow(tau,-1)/12.0+(c[0]*c[1]+2*b[0]*b[1])*etas*s;
Ak[0][2] = rho*s*pow(tau,-1)/12.0+(c[0]*c[2]+2*b[0]*b[2])*etas*s;
Ak[0][3] = b[0]*c[0]*etas*s;
Ak[0][4] = c[0]*b[1]*etas*s;
Ak[0][5] = c[0]*b[2]*etas*s;
Ak[0][6] = -b[0]*s/3.0;
Ak[0][7] = -b[0]*s/3.0;
Ak[0][8] = -b[0]*s/3.0;
Ak[0][9] = b[0]*s/3.0;
Ak[0][10] = b[0]*s/3.0;
・
・
・
Ak[20][17] = 0;
Ak[20][18] = s*(pow(tau,-1)*lambda+1)/12.0;
Ak[20][19] = s*(pow(tau,-1)*lambda+1)/12.0;
Ak[20][20] = s*(pow(tau,-1)*lambda+1)/6.0;
 
}
void make_fk(double *fk,double *uk,double *b,double tau,double s){
fk[0] = s*(uk[2]*rho*pow(tau,-1)+b[2])/12.0+s*(uk[1]*rho*pow(tau,-1)+b[1])/
    12.0+s*(uk[0]*rho*pow(tau,-1)+b[0])/6.0;
fk[1] = s*(uk[2]*rho*pow(tau,-1)+b[2])/12.0+s*(uk[1]*rho*pow(tau,-1)+b[1])/
    6.0+s*(uk[0]*rho*pow(tau,-1)+b[0])/12.0;
fk[2] = s*(uk[2]*rho*pow(tau,-1)+b[2])/6.0+s*(uk[1]*rho*pow(tau,-1)+b[1])/
    12.0+s*(uk[0]*rho*pow(tau,-1)+b[0])/12.0;
・
・
・
    6.0+s*(uk[18]*pow(tau,-1)*lambda+b[18])/12.0;
fk[20] = s*(uk[20]*pow(tau,-1)*lambda+b[20])/
    6.0+s*(uk[19]*pow(tau,-1)*lambda+b[19])/
    12.0+s*(uk[18]*pow(tau,-1)*lambda+b[18])/12.0;
 
}