Artifact 9d8d813cd9572a966be72ce20bb4709c0fd99f2d4674bfee00f098000fa2ba8a:
- File
r34/lib/fide.log
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 107670) [annotate] [blame] [check-ins using] [more...]
REDUCE 3.4, 15-Jul-91 ... 1: SF(1) := 1 SF(2) := 1 SF(3) := 1 (FIDE) %*********************************************************************** %***** ***** %***** Package F I D E --- T e s t Examples Ver. 1.1 05/1991 ***** %***** ***** %*********************************************************************** %*********************************************************************** %***** ***** %***** T e s t Examples --- Module E X P R E S ***** %***** ***** %*********************************************************************** let cos th**2=1 - sin th**2, cos fi**2=1 - sin fi**2; factor df; on rat; for all x,y let diff(x,y)=df(x,y); depend u,r,th,fi; depend v,r,th,fi; depend f,r,th,fi; depend w,r,th,fi; % Spherical coordinate system scalefactors 3,r*sin th*cos fi,r*sin th*sin fi,r*cos th,r,th,fi; SF(1) := 1 SF(2) := R SF(3) := SIN(TH)*R tensor a1,a2,a3,a4,a5; vectors u,v; dyads w; a1:=grad f; A1:= ( DF(F,R) , DF(F,TH) ---------- , R DF(F,FI) ----------- ) SIN(TH)*R a2:=div u; DF(U(3),FI) DF(U(2),TH) A2:=------------- + ------------- + DF(U(1),R) SIN(TH)*R R U(2)*COS(TH) + 2*U(1)*SIN(TH) + ------------------------------- SIN(TH)*R a3:=curl v; DF(V(3),TH) - DF(V(2),FI) V(3)*COS(TH) A3:= ( ------------- + ---------------- + -------------- , R SIN(TH)*R SIN(TH)*R DF(V(1),FI) - V(3) - DF(V(3),R) + ------------- + --------- , SIN(TH)*R R - DF(V(1),TH) V(2) DF(V(2),R) + ---------------- + ------ ) R R a4:=lapl v; - 2*DF(V(3),FI) - 2*DF(V(2),TH) A4:= ( ------------------ + ------------------ + DF(V(1),R,2) 2 2 SIN(TH)*R R 2*DF(V(1),R) DF(V(1),TH,2) DF(V(1),TH)*COS(TH) + -------------- + --------------- + --------------------- R 2 2 R SIN(TH)*R DF(V(1),FI,2) 2*(V(2)*COS(TH) + V(1)*SIN(TH)) + --------------- - --------------------------------- , 2 2 2 SIN(TH) *R SIN(TH)*R - 2*DF(V(3),FI)*COS(TH) 2*DF(V(2),R) -------------------------- + DF(V(2),R,2) + -------------- 2 2 R SIN(TH) *R DF(V(2),TH,2) DF(V(2),TH)*COS(TH) DF(V(2),FI,2) + --------------- + --------------------- + --------------- 2 2 2 2 R SIN(TH)*R SIN(TH) *R 2*DF(V(1),TH) - V(2) + --------------- + ------------- , 2 2 2 R SIN(TH) *R 2*DF(V(3),R) DF(V(3),TH,2) DF(V(3),R,2) + -------------- + --------------- R 2 R DF(V(3),TH)*COS(TH) DF(V(3),FI,2) + --------------------- + --------------- 2 2 2 SIN(TH)*R SIN(TH) *R 2*DF(V(2),FI)*COS(TH) 2*DF(V(1),FI) - V(3) + ----------------------- + --------------- + ------------- 2 2 2 2 2 SIN(TH) *R SIN(TH)*R SIN(TH) *R ) a3:=2*a3+a4; 2*DF(V(3),TH) - 2*DF(V(3),FI) - 2*DF(V(2),TH) A3:= ( --------------- + ------------------ + ------------------ R 2 2 SIN(TH)*R R - 2*DF(V(2),FI) 2*DF(V(1),R) + ------------------ + DF(V(1),R,2) + -------------- SIN(TH)*R R DF(V(1),TH,2) DF(V(1),TH)*COS(TH) DF(V(1),FI,2) + --------------- + --------------------- + --------------- 2 2 2 2 R SIN(TH)*R SIN(TH) *R 2*(V(3)*COS(TH)*R - V(2)*COS(TH) - V(1)*SIN(TH)) + -------------------------------------------------- , 2 SIN(TH)*R - 2*DF(V(3),FI)*COS(TH) - 2*DF(V(3),R) + -------------------------- + DF(V(2),R,2) 2 2 SIN(TH) *R 2*DF(V(2),R) DF(V(2),TH,2) DF(V(2),TH)*COS(TH) + -------------- + --------------- + --------------------- R 2 2 R SIN(TH)*R DF(V(2),FI,2) 2*DF(V(1),TH) 2*DF(V(1),FI) + --------------- + --------------- + --------------- 2 2 2 SIN(TH)*R SIN(TH) *R R 2 - 2*V(3)*SIN(TH) *R - V(2) + ----------------------------- , 2 2 SIN(TH) *R 2*DF(V(3),R) DF(V(3),TH,2) DF(V(3),R,2) + -------------- + --------------- R 2 R DF(V(3),TH)*COS(TH) DF(V(3),FI,2) + --------------------- + --------------- + 2*DF(V(2),R) 2 2 2 SIN(TH)*R SIN(TH) *R 2*DF(V(2),FI)*COS(TH) - 2*DF(V(1),TH) + ----------------------- + ------------------ 2 2 R SIN(TH) *R 2 2*DF(V(1),FI) - V(3) + 2*V(2)*SIN(TH) *R + --------------- + ----------------------------- ) 2 2 2 SIN(TH)*R SIN(TH) *R a5:=lapl f; 2*DF(F,R) DF(F,TH,2) DF(F,TH)*COS(TH) A5:=DF(F,R,2) + ----------- + ------------ + ------------------ R 2 2 R SIN(TH)*R DF(F,FI,2) + ------------- 2 2 SIN(TH) *R a1:=a1+div w; DF(W(3,1),FI) DF(W(2,1),TH) A1:= ( --------------- + --------------- + DF(W(1,1),R) + DF(F,R) + ( SIN(TH)*R R - W(3,3)*SIN(TH) - W(2,2)*SIN(TH) + W(2,1)*COS(TH) + 2*W(1,1)*SIN(TH))/(SIN(TH)*R) , DF(W(3,2),FI) DF(W(2,2),TH) DF(F,TH) --------------- + --------------- + DF(W(1,2),R) + ---------- SIN(TH)*R R R + ( - W(3,3)*COS(TH) + W(2,2)*COS(TH) + W(2,1)*SIN(TH) + 2*W(1,2)*SIN(TH))/(SIN(TH)*R) , DF(W(3,3),FI) DF(W(2,3),TH) DF(F,FI) --------------- + --------------- + DF(W(1,3),R) + ----------- SIN(TH)*R R SIN(TH)*R + (W(3,2)*COS(TH) + W(3,1)*SIN(TH) + W(2,3)*COS(TH) + 2*W(1,3)*SIN(TH))/(SIN(TH)*R) ) a1:=u.dyad((a,0,1),(1,b,3),(0,c,d)); A1:= ( U(2) + U(1)*A , U(3)*C + U(2)*B , U(3)*D + 3*U(2) + U(1) ) a2:=vect(a,b,c); A2:= ( A , B , C ) a1.a2; 2 2 U(3)*B*C + U(3)*C*D + U(2)*A + U(2)*B + 3*U(2)*C + U(1)*A + U(1)*C % Scalar product u.v; U(3)*V(3) + U(2)*V(2) + U(1)*V(1) % Vector product u?v; ( - U(3)*V(2) + U(2)*V(3) , U(3)*V(1) - U(1)*V(3) , - U(2)*V(1) + U(1)*V(2) ) % Dyadic u&v; ( ( U(1)*V(1) , U(1)*V(2) , U(1)*V(3) ) , ( U(2)*V(1) , U(2)*V(2) , U(2)*V(3) ) , ( U(3)*V(1) , U(3)*V(2) , U(3)*V(3) ) ) % Directional derivative dirdf(u,v); DF(V(1),TH)*U(2) DF(V(1),FI)*U(3) ( DF(V(1),R)*U(1) + ------------------ + ------------------ R SIN(TH)*R U(3)*V(3) + U(2)*V(2) - ----------------------- , R DF(V(2),TH)*U(2) DF(V(2),FI)*U(3) DF(V(2),R)*U(1) + ------------------ + ------------------ R SIN(TH)*R - U(3)*V(3)*COS(TH) + U(2)*V(1)*SIN(TH) + ------------------------------------------ , SIN(TH)*R DF(V(3),TH)*U(2) DF(V(3),FI)*U(3) DF(V(3),R)*U(1) + ------------------ + ------------------ R SIN(TH)*R U(3)*(V(2)*COS(TH) + V(1)*SIN(TH)) + ------------------------------------ ) SIN(TH)*R clear a1,a2,a3,a4,a5,u,v,w; for all x,y clear diff(x,y); clear cos th**2, cos fi**2; remfac df; off rat; scalefactors 3,x,y,z,x,y,z; SF(1) := 1 SF(2) := 1 SF(3) := 1 %*********************************************************************** %***** ***** %***** T e s t Examples --- Module I I M E T ***** %***** ***** %*********************************************************************** % Example I.1 - 1-D Lagrangian Hydrodynamics off exp; factor diff; on rat,eqfu; % Declare which indexes will be given to coordinates coordinates x,t into j,m; % Declares uniform grid in x coordinate grid uniform,x; % Declares dependencies of functions on coordinates dependence eta(t,x),v(t,x),eps(t,x),p(t,x); % Declares p as known function given p; same eta,v,p; iim a, eta,diff(eta,t)-eta*diff(v,x)=0, v,diff(v,t)+eta/ro*diff(p,x)=0, eps,diff(eps,t)+eta*p/ro*diff(v,x)=0; ***************************** ***** Program ***** IIMET Ver 1.1 17/05/91 ***************************** Partial Differential Equations ============================== - DIFF(V,X)*ETA + DIFF(ETA,T) = 0 DIFF(P,X)*ETA --------------- + DIFF(V,T) = 0 RO DIFF(V,X)*P*ETA ----------------- + DIFF(EPS,T) = 0 RO Backtracking needed in grid optimalization 0 interpolations are needed in X coordinate Equation for ETA variable is integrated in HALF grid point Equation for V variable is integrated in HALF grid point Equation for EPS variable is integrated in HALF grid point 0 interpolations are needed in T coordinate Equation for ETA variable is integrated in HALF grid point Equation for V variable is integrated in HALF grid point Equation for EPS variable is integrated in HALF grid point Equations after Discretization Using IIM : ========================================== ( - ((V(J + 1,M + 1) - V(J,M + 1))*(ETA(J + 1,M + 1) + ETA(J,M + 1)) + (V(J + 1,M) - V(J,M))*(ETA(J + 1,M) + ETA(J,M))) *(HT(M + 1) + HT(M)) + 4*(ETA(J + 1,M + 1) - ETA(J + 1,M) + ETA(J,M + 1) - ETA(J,M))*HX) /(4*(HT(M + 1) + HT(M))*HX) = 0 (((P(J + 1,M + 1) - P(J,M + 1))*(ETA(J + 1,M + 1) + ETA(J,M + 1)) + (P(J + 1,M) - P(J,M))*(ETA(J + 1,M) + ETA(J,M))) *(HT(M + 1) + HT(M)) + 4*(V(J + 1,M + 1) - V(J + 1,M) + V(J,M + 1) - V(J,M))*RO*HX)/(4 *(HT(M + 1) + HT(M))*RO*HX) = 0 (((P(J + 1,M + 1)*ETA(J + 1,M + 1) + P(J,M + 1)*ETA(J,M + 1)) *(V(J + 1,M + 1) - V(J,M + 1)) + (P(J + 1,M)*ETA(J + 1,M) + P(J,M)*ETA(J,M))*(V(J + 1,M) - V(J,M))) *(HT(M + 1) + HT(M)) + 4 *(EPS(J + 1,M + 1) - EPS(J + 1,M) + EPS(J,M + 1) - EPS(J,M))*RO*HX)/ (4*(HT(M + 1) + HT(M))*RO*HX) = 0 clear a; clearsame; cleargiven; %*********************************************************************** % Example I.2 - How other functions (here sin, cos) can be used in % discretized terms diffunc sin,cos; difmatch all,diff(u*sin x,x),u=one,2,(u(i+1)*sin x(i+1)-u(i-1) *sin x(i-1))/(dim1+dip1), u=half,0,(u(i+1/2)*sin x(i+1/2)-u(i-1/2)*sin x(i-1/2)) /di; difmatch all,cos x*diff(u,x,2),u=one,0,cos x i*(u(i+1)-2*u(i)+u(i-1)) /di^2, u=half,3,(u(i+3/2)-u(i+1/2))/dip2/2 - (u(i-1/2)-u(i-3/2))/dim2/2; off exp; coordinates x,t into j,m; grid uniform,x,t; dependence u(x,t),v(x,t); iim a,u,diff(u,t)+diff(u,x)+cos x*diff(v,x,2)=0, v,diff(v,t)+diff(sin x*u,x)=0; ***************************** ***** Program ***** IIMET Ver 1.1 17/05/91 ***************************** Partial Differential Equations ============================== DIFF(U,T) + DIFF(U,X) + DIFF(V,X,2)*COS(X) = 0 DIFF(SIN(X)*U,X) + DIFF(V,T) = 0 0 interpolations are needed in X coordinate Equation for U variable is integrated in HALF grid point Equation for V variable is integrated in HALF grid point 0 interpolations are needed in T coordinate Equation for U variable is integrated in HALF grid point Equation for V variable is integrated in HALF grid point Equations after Discretization Using IIM : ========================================== (((U(J + 1,M + 1) + U(J + 1,M))*HT + (U(J + 1,M + 1) - U(J + 1,M) + U(J,M + 1) - U(J,M))*HX 2*J - 1 - (U(J,M + 1) + U(J,M))*HT)*HX + (V(---------,M + 1) 2 2*J - 1 2*J + 3 2*J + 3 + V(---------,M) + V(---------,M + 1) + V(---------,M) 2 2 2 2*J + 1 2*J + 1 2*J + 1 - 2*V(---------,M + 1) - 2*V(---------,M))*COS(X(---------))*HT) 2 2 2 2 /(2*HX *HT) = 0 (((U(J + 1,M + 1) + U(J + 1,M))*SIN(X(J + 1)) - (U(J,M + 1) + U(J,M))*SIN(X(J)))*HT 2*J + 1 2*J + 1 + 2*(V(---------,M + 1) - V(---------,M))*HX)/(2*HX*HT) = 0 2 2 clear a; %*********************************************************************** % Example I.3 - Schrodinger equation factor diff; coordinates t,x into m,j; grid uniform,x,t; dependence ur(x,t),ui(x,t); same ui,ur; iim a,ur,-diff(ui,t)+1/2*diff(ur,x,2)+(ur**2+ui**2)*ur=0, ui,diff(ur,t)+1/2*diff(ui,x,2)+(ur**2+ui**2)*ui=0; ***************************** ***** Program ***** IIMET Ver 1.1 17/05/91 ***************************** Partial Differential Equations ============================== DIFF(UR,X,2) 2 2 -------------- - DIFF(UI,T) = - UR*(UR + UI ) 2 DIFF(UI,X,2) 2 2 DIFF(UR,T) + -------------- = - UI*(UR + UI ) 2 0 interpolations are needed in T coordinate Equation for UR variable is integrated in HALF grid point Equation for UI variable is integrated in HALF grid point 0 interpolations are needed in X coordinate Equation for UR variable is integrated in ONE grid point Equation for UI variable is integrated in ONE grid point Equations after Discretization Using IIM : ========================================== ((UR(M + 1,J - 1) + UR(M + 1,J + 1) - 2*UR(M + 1,J) + UR(M,J - 1) 2 2 + UR(M,J + 1) - 2*UR(M,J))*HT - 4*UI(M + 1,J)*HX + 4*UI(M,J)*HX ) 2 3 2 /(4*HX *HT) = - (UR(M + 1,J) + UR(M + 1,J)*UI(M + 1,J) 3 2 + UR(M,J) + UR(M,J)*UI(M,J) )/2 ((UI(M + 1,J - 1) + UI(M + 1,J + 1) - 2*UI(M + 1,J) + UI(M,J - 1) 2 2 + UI(M,J + 1) - 2*UI(M,J))*HT + 4*UR(M + 1,J)*HX - 4*UR(M,J)*HX ) 2 2 2 /(4*HX *HT) = - (UR(M + 1,J) *UI(M + 1,J) + UR(M,J) *UI(M,J) 3 3 + UI(M + 1,J) + UI(M,J) )/2 clear a; clearsame; %*********************************************************************** % Example I.4 - Vector calculus in p.d.e. input % cooperation with expres module % 2-D hydrodynamics scalefactors 2,x,y,x,y; SF(1) := 1 SF(2) := 1 vectors u; off exp,twogrid; on eqfu; factor diff,ht,hx,hy; coordinates x,y,t into j,i,m; grid uniform,x,y,t; dependence n(t,x,y),u(t,x,y),p(t,x,y); iim a,n,diff(n,t)+u.grad n+n*div u=0, u,m*n*(diff(u,t)+u.grad u)+grad p=vect(0,0), p,3/2*(diff(p,t)+u.grad p)+5/2*p*div u=0; ***************************** ***** Program ***** IIMET Ver 1.1 17/05/91 ***************************** Partial Differential Equations ============================== DIFF(N,T) + DIFF(N,X)*U1 + DIFF(N,Y)*U2 + DIFF(U1,X)*N + DIFF(U2,Y)*N = 0 DIFF(P,X) + DIFF(U1,T)*M*N + DIFF(U1,X)*M*N*U1 + DIFF(U1,Y)*M*N*U2 = 0 DIFF(P,Y) + DIFF(U2,T)*M*N + DIFF(U2,X)*M*N*U1 + DIFF(U2,Y)*M*N*U2 = 0 3*DIFF(P,T) 3*DIFF(P,X)*U1 3*DIFF(P,Y)*U2 ------------- + ---------------- + ---------------- 2 2 2 5*DIFF(U1,X)*P 5*DIFF(U2,Y)*P + ---------------- + ---------------- = 0 2 2 0 interpolations are needed in X coordinate Equation for N variable is integrated in HALF grid point Equation for U1 variable is integrated in HALF grid point Equation for U2 variable is integrated in HALF grid point Equation for P variable is integrated in HALF grid point 0 interpolations are needed in Y coordinate Equation for N variable is integrated in HALF grid point Equation for U1 variable is integrated in HALF grid point Equation for U2 variable is integrated in HALF grid point Equation for P variable is integrated in HALF grid point 0 interpolations are needed in T coordinate Equation for N variable is integrated in HALF grid point Equation for U1 variable is integrated in HALF grid point Equation for U2 variable is integrated in HALF grid point Equation for P variable is integrated in HALF grid point Equations after Discretization Using IIM : ========================================== -1 -1 (HY *HX *((U1(J + 1,I + 1,M + 1)*HY + U2(J + 1,I + 1,M + 1)*HX) *N(J + 1,I + 1,M + 1) + (U1(J + 1,I + 1,M)*HY + U2(J + 1,I + 1,M)*HX) *N(J + 1,I + 1,M) + (U1(J + 1,I,M + 1)*HY - U2(J + 1,I,M + 1)*HX) *N(J + 1,I,M + 1) + (U1(J + 1,I,M)*HY - U2(J + 1,I,M)*HX)*N(J + 1,I,M) - (U1(J,I + 1,M + 1)*HY - U2(J,I + 1,M + 1)*HX) *N(J,I + 1,M + 1) - (U1(J,I + 1,M)*HY - U2(J,I + 1,M)*HX)*N(J,I + 1,M) - (U1(J,I,M + 1)*HY + U2(J,I,M + 1)*HX)*N(J,I,M + 1) -1 - (U1(J,I,M)*HY + U2(J,I,M)*HX)*N(J,I,M)))/4 + (HT *( N(J + 1,I + 1,M + 1) - N(J + 1,I + 1,M) + N(J + 1,I,M + 1) - N(J + 1,I,M) + N(J,I + 1,M + 1) - N(J,I + 1,M) + N(J,I,M + 1) - N(J,I,M)))/4 = 0 -1 (HX *((P(J + 1,I + 1,M + 1) + P(J + 1,I + 1,M) + P(J + 1,I,M + 1) + P(J + 1,I,M)) - (P(J,I + 1,M + 1) + P(J,I + 1,M) + P(J,I,M + 1) + P(J,I,M)) -1 -1 ))/4 + (HY *HX *((( N(J + 1,I + 1,M + 1)*U1(J + 1,I + 1,M + 1) + N(J,I + 1,M + 1)*U1(J,I + 1,M + 1)) *(U1(J + 1,I + 1,M + 1) - U1(J,I + 1,M + 1)) + ( N(J + 1,I + 1,M)*U1(J + 1,I + 1,M) + N(J,I + 1,M)*U1(J,I + 1,M)) *(U1(J + 1,I + 1,M) - U1(J,I + 1,M)) + ( N(J + 1,I,M + 1)*U1(J + 1,I,M + 1) + N(J,I,M + 1)*U1(J,I,M + 1)) *(U1(J + 1,I,M + 1) - U1(J,I,M + 1)) + ( N(J + 1,I,M)*U1(J + 1,I,M) + N(J,I,M)*U1(J,I,M)) *(U1(J + 1,I,M) - U1(J,I,M)))*HY + (( N(J + 1,I + 1,M + 1)*U2(J + 1,I + 1,M + 1) + N(J + 1,I,M + 1)*U2(J + 1,I,M + 1)) *(U1(J + 1,I + 1,M + 1) - U1(J + 1,I,M + 1)) + ( N(J + 1,I + 1,M)*U2(J + 1,I + 1,M) + N(J + 1,I,M)*U2(J + 1,I,M)) *(U1(J + 1,I + 1,M) - U1(J + 1,I,M)) + ( N(J,I + 1,M + 1)*U2(J,I + 1,M + 1) + N(J,I,M + 1)*U2(J,I,M + 1)) *(U1(J,I + 1,M + 1) - U1(J,I,M + 1)) + ( N(J,I + 1,M)*U2(J,I + 1,M) + N(J,I,M)*U2(J,I,M)) -1 *(U1(J,I + 1,M) - U1(J,I,M)))*HX)*M)/8 + (HT *( (N(J + 1,I + 1,M + 1) + N(J + 1,I + 1,M)) *(U1(J + 1,I + 1,M + 1) - U1(J + 1,I + 1,M)) + (N(J + 1,I,M + 1) + N(J + 1,I,M)) *(U1(J + 1,I,M + 1) - U1(J + 1,I,M)) + (N(J,I + 1,M + 1) + N(J,I + 1,M)) *(U1(J,I + 1,M + 1) - U1(J,I + 1,M)) + (N(J,I,M + 1) + N(J,I,M))*(U1(J,I,M + 1) - U1(J,I,M)))*M)/8 = 0 -1 (HY *((N(J + 1,I + 1,M + 1)*U2(J + 1,I + 1,M + 1) + N(J + 1,I,M + 1)*U2(J + 1,I,M + 1)) *(U2(J + 1,I + 1,M + 1) - U2(J + 1,I,M + 1)) + ( N(J + 1,I + 1,M)*U2(J + 1,I + 1,M) + N(J + 1,I,M)*U2(J + 1,I,M)) *(U2(J + 1,I + 1,M) - U2(J + 1,I,M)) + ( N(J,I + 1,M + 1)*U2(J,I + 1,M + 1) + N(J,I,M + 1)*U2(J,I,M + 1)) *(U2(J,I + 1,M + 1) - U2(J,I,M + 1)) + (N(J,I + 1,M)*U2(J,I + 1,M) + N(J,I,M)*U2(J,I,M)) -1 *(U2(J,I + 1,M) - U2(J,I,M)))*M)/8 + (HX *(( N(J + 1,I + 1,M + 1)*U1(J + 1,I + 1,M + 1) + N(J,I + 1,M + 1)*U1(J,I + 1,M + 1)) *(U2(J + 1,I + 1,M + 1) - U2(J,I + 1,M + 1)) + ( N(J + 1,I + 1,M)*U1(J + 1,I + 1,M) + N(J,I + 1,M)*U1(J,I + 1,M)) *(U2(J + 1,I + 1,M) - U2(J,I + 1,M)) + ( N(J + 1,I,M + 1)*U1(J + 1,I,M + 1) + N(J,I,M + 1)*U1(J,I,M + 1)) *(U2(J + 1,I,M + 1) - U2(J,I,M + 1)) + (N(J + 1,I,M)*U1(J + 1,I,M) + N(J,I,M)*U1(J,I,M)) -1 -1 *(U2(J + 1,I,M) - U2(J,I,M)))*M)/8 + (HY *HT *(( (N(J + 1,I + 1,M + 1) + N(J + 1,I + 1,M)) *(U2(J + 1,I + 1,M + 1) - U2(J + 1,I + 1,M)) + (N(J + 1,I,M + 1) + N(J + 1,I,M)) *(U2(J + 1,I,M + 1) - U2(J + 1,I,M)) + (N(J,I + 1,M + 1) + N(J,I + 1,M)) *(U2(J,I + 1,M + 1) - U2(J,I + 1,M)) + (N(J,I,M + 1) + N(J,I,M))*(U2(J,I,M + 1) - U2(J,I,M)))*M *HY + 2*(P(J + 1,I + 1,M + 1) + P(J + 1,I + 1,M))*HT - 2*(P(J + 1,I,M + 1) + P(J + 1,I,M))*HT + 2*(P(J,I + 1,M + 1) + P(J,I + 1,M))*HT - 2*(P(J,I,M + 1) + P(J,I,M))*HT))/8 = 0 -1 -1 (HY *HX *(5*((P(J + 1,I + 1,M + 1) + P(J + 1,I,M + 1)) *(U2(J + 1,I + 1,M + 1) - U2(J + 1,I,M + 1)) + (P(J + 1,I + 1,M) + P(J + 1,I,M)) *(U2(J + 1,I + 1,M) - U2(J + 1,I,M)) + (P(J,I + 1,M + 1) + P(J,I,M + 1)) *(U2(J,I + 1,M + 1) - U2(J,I,M + 1)) + (P(J,I + 1,M) + P(J,I,M))*(U2(J,I + 1,M) - U2(J,I,M))) *HX + 3*((P(J + 1,I + 1,M + 1) - P(J + 1,I,M + 1)) *(U2(J + 1,I + 1,M + 1) + U2(J + 1,I,M + 1)) + (P(J + 1,I + 1,M) - P(J + 1,I,M)) *(U2(J + 1,I + 1,M) + U2(J + 1,I,M)) + (P(J,I + 1,M + 1) - P(J,I,M + 1)) *(U2(J,I + 1,M + 1) + U2(J,I,M + 1)) + (P(J,I + 1,M) - P(J,I,M))*(U2(J,I + 1,M) + U2(J,I,M))) *HX + 5*((P(J + 1,I + 1,M + 1) + P(J,I + 1,M + 1)) *(U1(J + 1,I + 1,M + 1) - U1(J,I + 1,M + 1)) + (P(J + 1,I + 1,M) + P(J,I + 1,M)) *(U1(J + 1,I + 1,M) - U1(J,I + 1,M)) + (P(J + 1,I,M + 1) + P(J,I,M + 1)) *(U1(J + 1,I,M + 1) - U1(J,I,M + 1)) + (P(J + 1,I,M) + P(J,I,M))*(U1(J + 1,I,M) - U1(J,I,M))) *HY + 3*((P(J + 1,I + 1,M + 1) - P(J,I + 1,M + 1)) *(U1(J + 1,I + 1,M + 1) + U1(J,I + 1,M + 1)) + (P(J + 1,I + 1,M) - P(J,I + 1,M)) *(U1(J + 1,I + 1,M) + U1(J,I + 1,M)) + (P(J + 1,I,M + 1) - P(J,I,M + 1)) *(U1(J + 1,I,M + 1) + U1(J,I,M + 1)) + (P(J + 1,I,M) - P(J,I,M))*(U1(J + 1,I,M) + U1(J,I,M))) -1 *HY))/16 + (3*HT *(P(J + 1,I + 1,M + 1) - P(J + 1,I + 1,M) + P(J + 1,I,M + 1) - P(J + 1,I,M) + P(J,I + 1,M + 1) - P(J,I + 1,M) + P(J,I,M + 1) - P(J,I,M)))/ 8 = 0 clear a,u; %*********************************************************************** % Example I.5 - 1-D hydrodynamics up to 3-rd moments (heat flow) coordinates x,t into j,m; grid uniform,x,t; dependence n(x,t),u(x,t),tt(x,t),p(x,t),q(x,t); iim a, n,diff(n,t)+u*diff(n,x)+diff(u,x)=0, u,n*m*(diff(u,t)+u*diff(u,x))+k*diff(n*tt,x)+diff(p,x)=0, tt,3/2*k*n*(diff(tt,t)+u*diff(tt,x))+n*k*tt*diff(u,x)+1/2*p *diff(u,x)+diff(q,x)=0, p,diff(p,t)+u*diff(p,x)+p*diff(u,x)+n*k*tt*diff(u,x)+2/5*diff(q,x) =0, q,diff(q,t)+u*diff(q,x)+q*diff(u,x)+5/2*n*k**2*tt/m*diff(tt,x)+n*k *tt*diff(p,x)-p*diff(p,x)=0; ***************************** ***** Program ***** IIMET Ver 1.1 17/05/91 ***************************** Partial Differential Equations ============================== DIFF(N,T) + DIFF(N,X)*U + DIFF(U,X) = 0 DIFF(N*TT,X)*K + DIFF(P,X) + DIFF(U,T)*M*N + DIFF(U,X)*M*N*U = 0 DIFF(U,X)*(2*K*N*TT + P) 3*DIFF(TT,T)*K*N DIFF(Q,X) + -------------------------- + ------------------ 2 2 3*DIFF(TT,X)*K*N*U + -------------------- = 0 2 2*DIFF(Q,X) DIFF(P,T) + DIFF(P,X)*U + ------------- + DIFF(U,X)*(K*N*TT + P) 5 = 0 DIFF(P,X)*(K*N*TT - P) + DIFF(Q,T) + DIFF(Q,X)*U + DIFF(U,X)*Q 2 5*DIFF(TT,X)*K *N*TT + ---------------------- = 0 2*M 0 interpolations are needed in X coordinate Equation for N variable is integrated in HALF grid point Equation for U variable is integrated in HALF grid point Equation for TT variable is integrated in HALF grid point Equation for P variable is integrated in HALF grid point Equation for Q variable is integrated in HALF grid point 0 interpolations are needed in T coordinate Equation for N variable is integrated in HALF grid point Equation for U variable is integrated in HALF grid point Equation for TT variable is integrated in HALF grid point Equation for P variable is integrated in HALF grid point Equation for Q variable is integrated in HALF grid point Equations after Discretization Using IIM : ========================================== -1 (HX *(((N(J + 1,M + 1) - N(J,M + 1))*(U(J + 1,M + 1) + U(J,M + 1)) + (N(J + 1,M) - N(J,M))*(U(J + 1,M) + U(J,M))) + 2*(U(J + 1,M + 1) + U(J + 1,M)) - 2*(U(J,M + 1) + U(J,M)))) -1 HT *(N(J + 1,M + 1) - N(J + 1,M) + N(J,M + 1) - N(J,M)) /4 + ---------------------------------------------------------- 2 = 0 -1 (HX *(((N(J + 1,M + 1)*U(J + 1,M + 1) + N(J,M + 1)*U(J,M + 1)) *(U(J + 1,M + 1) - U(J,M + 1)) + (N(J + 1,M)*U(J + 1,M) + N(J,M)*U(J,M))*(U(J + 1,M) - U(J,M)) )*M + 2*(N(J + 1,M + 1)*TT(J + 1,M + 1) + N(J + 1,M)*TT(J + 1,M) - N(J,M + 1)*TT(J,M + 1) - N(J,M)*TT(J,M))*K + 2*(P(J + 1,M + 1) + P(J + 1,M)) - 2*(P(J,M + 1) + P(J,M)))) -1 /4 + (HT *((N(J + 1,M + 1) + N(J + 1,M)) *(U(J + 1,M + 1) - U(J + 1,M)) + (N(J,M + 1) + N(J,M))*(U(J,M + 1) - U(J,M)))*M)/4 = 0 -1 (HX *(3*((N(J + 1,M + 1)*U(J + 1,M + 1) + N(J,M + 1)*U(J,M + 1)) *(TT(J + 1,M + 1) - TT(J,M + 1)) + (N(J + 1,M)*U(J + 1,M) + N(J,M)*U(J,M)) *(TT(J + 1,M) - TT(J,M)))*K + 2*( (N(J + 1,M + 1)*TT(J + 1,M + 1) + N(J,M + 1)*TT(J,M + 1)) *(U(J + 1,M + 1) - U(J,M + 1)) + (N(J + 1,M)*TT(J + 1,M) + N(J,M)*TT(J,M)) *(U(J + 1,M) - U(J,M)))*K + ( (P(J + 1,M + 1) + P(J,M + 1))*(U(J + 1,M + 1) - U(J,M + 1)) + (P(J + 1,M) + P(J,M))*(U(J + 1,M) - U(J,M))) + 4*(Q(J + 1,M + 1) + Q(J + 1,M)) - 4*(Q(J,M + 1) + Q(J,M)))) -1 /8 + (3*HT *((N(J + 1,M + 1) + N(J + 1,M)) *(TT(J + 1,M + 1) - TT(J + 1,M)) + (N(J,M + 1) + N(J,M))*(TT(J,M + 1) - TT(J,M)))*K)/8 = 0 -1 (HX *(5*(((N(J + 1,M + 1)*TT(J + 1,M + 1) + N(J,M + 1)*TT(J,M + 1)) *(U(J + 1,M + 1) - U(J,M + 1)) + (N(J + 1,M)*TT(J + 1,M) + N(J,M)*TT(J,M)) *(U(J + 1,M) - U(J,M)))*K + (P(J + 1,M + 1) + P(J,M + 1)) *(U(J + 1,M + 1) - U(J,M + 1)) + (P(J + 1,M + 1) - P(J,M + 1))*(U(J + 1,M + 1) + U(J,M + 1)) + (P(J + 1,M) + P(J,M))*(U(J + 1,M) - U(J,M)) + (P(J + 1,M) - P(J,M))*(U(J + 1,M) + U(J,M))) + 4*(Q(J + 1,M + 1) + Q(J + 1,M)) - 4*(Q(J,M + 1) + Q(J,M)))) -1 HT *(P(J + 1,M + 1) - P(J + 1,M) + P(J,M + 1) - P(J,M)) /20 + ---------------------------------------------------------- 2 = 0 -1 (HX *(2*((N(J + 1,M + 1)*TT(J + 1,M + 1) + N(J,M + 1)*TT(J,M + 1)) *(P(J + 1,M + 1) - P(J,M + 1)) + (N(J + 1,M)*TT(J + 1,M) + N(J,M)*TT(J,M)) *(P(J + 1,M) - P(J,M)))*K*M + 5*( (N(J + 1,M + 1)*TT(J + 1,M + 1) + N(J,M + 1)*TT(J,M + 1)) *(TT(J + 1,M + 1) - TT(J,M + 1)) + (N(J + 1,M)*TT(J + 1,M) + N(J,M)*TT(J,M)) 2 *(TT(J + 1,M) - TT(J,M)))*K - 2*( (P(J + 1,M + 1) + P(J,M + 1))*(P(J + 1,M + 1) - P(J,M + 1)) + (P(J + 1,M) + P(J,M))*(P(J + 1,M) - P(J,M)))*M + 2*( (Q(J + 1,M + 1) + Q(J,M + 1))*(U(J + 1,M + 1) - U(J,M + 1)) + (Q(J + 1,M) + Q(J,M))*(U(J + 1,M) - U(J,M)))*M + 2*( (Q(J + 1,M + 1) - Q(J,M + 1))*(U(J + 1,M + 1) + U(J,M + 1)) + (Q(J + 1,M) - Q(J,M))*(U(J + 1,M) + U(J,M)))*M))/(8*M) -1 HT *(Q(J + 1,M + 1) - Q(J + 1,M) + Q(J,M + 1) - Q(J,M)) + ---------------------------------------------------------- = 0 2 clear a; remfac diff,ht,hx,hy; on exp; off rat; %*********************************************************************** %***** ***** %***** T e s t Examples --- Module A P P R O X ***** %***** ***** %*********************************************************************** % Example A.1 coordinates x,t into j,n; maxorder t=2,x=3; functions u,v; approx( (u(n+1/2)-u(n-1/2))/ht=(v(n+1/2,j+1/2)-v(n+1/2,j-1/2) +v(n-1/2,j+1/2)-v(n-1/2,j-1/2))/(2*hx) ); Difference scheme approximates differential equation DF(U,T)=DF(V,X) with orders of approximation: 2 HX 2 HT % Example A.2 maxorder t=3,x=3; approx( (u(n+1)-u(n))/ht=(u(n+1,j+1/2)-u(n+1,j-1/2) +u(n,j+1/2)-u(n,j-1/2))/(2*hx) ); Difference scheme approximates differential equation DF(U,T)=DF(U,X) with orders of approximation: 2 HX HT % Example A.3 maxorder t=2,x=3; center t=1/2; approx( (u(n+1)-u(n))/ht=(v(n+1,j+1/2)-v(n+1,j-1/2) +v(n,j+1/2)-v(n,j-1/2))/(2*hx) ); Difference scheme approximates differential equation DF(U,T)=DF(V,X) with orders of approximation: 2 HX 2 HT % Example A.4 approx( u(n+1)/ht=(v(n+1,j+1/2)-v(n+1,j-1/2) +v(n,j+1/2)-v(n,j-1/2))/(2*hx) ); Reformulate difference scheme, grid steps remain in denominators Difference scheme approximates differential equation 0=DF(V,X) with orders of approximation: 2 HX -1 HT % Example A.5 maxorder t=3,x=3; approx( (u(n+1)-u(n))/ht=(u(n+1,j+1/2)-u(n+1,j-1/2))/hx); Difference scheme approximates differential equation DF(U,T)=DF(U,X) with orders of approximation: 2 HX HT % Example A.6 approx( (u(n+1)-u(n))/ht=(u(n+1/2,j+1/2)-u(n+1/2,j-1/2))/hx); Difference scheme approximates differential equation DF(U,T)=DF(U,X) with orders of approximation: 2 HX 2 HT % Example A.7; approx((u(n+1)-u(n))/ht=(u(n+1/2,j+1)-2*u(n+1/2,j)+u(n+1/2,j-1))/hx**2); Difference scheme approximates differential equation DF(U,T) =DF(U,X,2) with orders of approximation: 1 2 HT %*********************************************************************** %***** ***** %***** T e s t Examples --- Module C H A R P O L ***** %***** ***** %*********************************************************************** % Example C.1 coordinates t,x into i,j; grid uniform,t,x; let cos ax**2=1-sin ax**2; unfunc u,v; matrix aa(1,2),bb(2,2); aa(1,1):=(u(i+1)-u(i))/ht+(v(j+1)-v(j))/hx$ aa(1,2):=(v(i+1)-v(i))/ht+(u(j+1/2)-u(j-1/2))/hx$ bb:=ampmat aa; KX*HX AX := ------- 2 [HX 0 ] H1 := [ ] [0 HX] 2 2 H0 := MAT((HX,HT*(SIN(AX) - 2*SIN(AX)*COS(AX)*I - COS(AX) + 1)), ( - 2*SIN(AX)*I*HT,HX)) [ 2*SIN(AX)*HT*(SIN(AX) - COS(AX)*I) ] [ 1 ------------------------------------] [ HX ] BB := [ ] [ - 2*SIN(AX)*I*HT ] [------------------- 1 ] [ HX ] bb:=denotemat bb; [ 1 I*AI12 + AR12] BB := [ ] [I*AI21 1 ] factor lam; pol:=charpol bb; 2 POL := LAM - 2*LAM - I*AR12*AI21 + AI12*AI21 + 1 prdenot; 2 2*SIN(AX) *HT AR12 := --------------- HX - 2*SIN(AX)*COS(AX)*HT AI12 := ------------------------- HX - 2*SIN(AX)*HT AI21 := ----------------- HX cleardenot; clear aa,bb,pol; %*********************************************************************** % Example C.2 : Reprint Vorozcov, Ganza, Mazurik: Simvolno-cislennyj % interfejs. v zadacach ..., Novosibirsk 1986, p.47. unfunc u; matrix aa(1,1),bb(1,1); aa(1,1):=(u(i+1)-u(i))/ht+a*(u(j)-u(j-1))/hx$ bb:=ampmat aa; AX := KX*HX H1 := [HX] H0 := [ - SIN(AX)*A*I*HT + COS(AX)*A*HT - A*HT + HX] [ - SIN(AX)*A*I*HT + COS(AX)*A*HT - A*HT + HX ] BB := [----------------------------------------------] [ HX ] bb:=denotemat bb; BB := [I*AI11 + AR11] pol:=charpol bb; POL := LAM - (I*AI11 + AR11) prdenot; COS(AX)*A*HT - A*HT + HX AR11 := -------------------------- HX - SIN(AX)*A*HT AI11 := ----------------- HX cleardenot; clear aa,bb,pol; %*********************************************************************** % Example C.3 : Reprint Vorozcov, Ganza, Mazurik: Simvolno-cislennyj % interfejs. v zadacach ..., Novosibirsk 1986, p.52. coordinates t,x into m,j; unfunc u,r; matrix aa(1,2),bb(2,2); aa(1,1):=(r(m+1)-r(m))/ht+u0*(r(m+1,j+1)-r(m+1,j-1))/2/hx +r0*(u(m+1,j+1)-u(m+1,j-1))/2/hx$ aa(1,2):=(u(m+1)-u(m))/ht+u0*(u(m+1,j+1)-u(m+1,j-1))/2/hx +c0**2/r0*(r(m,j+1)-u(m,j-1))/2/hx$ bb:=ampmat aa; AX := KX*HX [ SIN(AX)*I*HT*R0 SIN(AX)*I*HT*U0 + HX] H1 := [ ] [2*R0*(SIN(AX)*I*HT*U0 + HX) 0 ] H0 := MAT((0,HX), 2 2 ( - SIN(AX)*I*HT*C0 + COS(AX)*HT*C0 + 2*HX*R0, 2 - (HT*C0 )*(SIN(AX)*I + COS(AX)))) 2 2 2 2 2 BB := MAT((( - SIN(AX) *HT *U0*C0 - SIN(AX)*COS(AX)*I*HT *U0*C0 2 - 2*SIN(AX)*I*HX*HT*U0*R0 - SIN(AX)*I*HX*HT*C0 2 2 + COS(AX)*HX*HT*C0 + 2*HX *R0)/(2*R0 2 2 2 2 2 2 *(SIN(AX) *HT *U0 + HX )),(HT*C0 *( - SIN(AX) *HT*U0 + SIN(AX)*COS(AX)*I*HT*U0 - SIN(AX)*I*HX 2 2 2 2 - COS(AX)*HX))/(2*R0*(SIN(AX) *HT *U0 + HX ))), 3 3 2 2 ((SIN(AX)*HT*(SIN(AX) *HT *U0 *C0 2 3 2 2 + SIN(AX) *COS(AX)*I*HT *U0 *C0 2 2 2 + 2*SIN(AX) *I*HX*HT *U0 *R0 2 2 2 + 2*SIN(AX) *I*HX*HT *U0*C0 2 2 - 2*SIN(AX)*COS(AX)*HX*HT *U0*C0 2 2 2 - 4*SIN(AX)*HX *HT*U0*R0 - SIN(AX)*HX *HT*C0 2 2 3 - COS(AX)*I*HX *HT*C0 - 2*I*HX *R0))/(2 4 4 4 2 2 2 2 4 *(SIN(AX) *HT *U0 + 2*SIN(AX) *HX *HT *U0 + HX )),( 4 4 2 2 3 4 2 2 SIN(AX) *HT *U0 *C0 - SIN(AX) *COS(AX)*I*HT *U0 *C0 3 3 3 3 3 2 - 2*SIN(AX) *I*HX*HT *U0 + 2*SIN(AX) *I*HX*HT *U0*C0 2 3 2 + 2*SIN(AX) *COS(AX)*HX*HT *U0*C0 2 2 2 2 2 2 2 2 + 2*SIN(AX) *HX *HT *U0 - SIN(AX) *HX *HT *C0 2 2 2 + SIN(AX)*COS(AX)*I*HX *HT *C0 3 4 - 2*SIN(AX)*I*HX *HT*U0 + 2*HX )/(2 4 4 4 2 2 2 2 4 *(SIN(AX) *HT *U0 + 2*SIN(AX) *HX *HT *U0 + HX )))) bb:=denotemat bb; [I*AI11 + AR11 I*AI12 + AR12] BB := [ ] [I*AI21 + AR21 I*AI22 + AR22] pol:=charpol bb; 2 POL := LAM - LAM*(I*AI11 + I*AI22 + AR11 + AR22) - I*AR12*AI21 - I*AI12*AR21 + I*AR11*AI22 + I*AI11*AR22 - AR12*AR21 + AI12*AI21 + AR11*AR22 - AI11*AI22 prdenot; 2 2 2 2 2 - SIN(AX) *HT *U0*C0 + COS(AX)*HX*HT*C0 + 2*HX *R0 AR11 := ------------------------------------------------------- 2 2 2 2 2*R0*(SIN(AX) *HT *U0 + HX ) 2 2 SIN(AX)*HT*( - COS(AX)*HT*U0*C0 - 2*HX*U0*R0 - HX*C0 ) AI11 := --------------------------------------------------------- 2 2 2 2 2*R0*(SIN(AX) *HT *U0 + HX ) 2 2 - (HT*C0 )*(SIN(AX) *HT*U0 + COS(AX)*HX) AR12 := ------------------------------------------- 2 2 2 2 2*R0*(SIN(AX) *HT *U0 + HX ) 2 SIN(AX)*HT*C0 *(COS(AX)*HT*U0 - HX) AI12 := ------------------------------------- 2 2 2 2 2*R0*(SIN(AX) *HT *U0 + HX ) 2 2 2 2 2 2 2 AR21 := (SIN(AX) *HT *(SIN(AX) *HT *U0 *C0 - 2*COS(AX)*HX*HT*U0*C0 2 2 2 - 4*HX *U0*R0 - HX *C0 ))/(2 4 4 4 2 2 2 2 4 *(SIN(AX) *HT *U0 + 2*SIN(AX) *HX *HT *U0 + HX )) 2 3 2 2 AI21 := (SIN(AX)*HT*(SIN(AX) *COS(AX)*HT *U0 *C0 2 2 2 2 2 2 + 2*SIN(AX) *HX*HT *U0 *R0 + 2*SIN(AX) *HX*HT *U0*C0 2 2 3 - COS(AX)*HX *HT*C0 - 2*HX *R0))/(2 4 4 4 2 2 2 2 4 *(SIN(AX) *HT *U0 + 2*SIN(AX) *HX *HT *U0 + HX )) 4 4 2 2 2 3 2 AR22 := (SIN(AX) *HT *U0 *C0 + 2*SIN(AX) *COS(AX)*HX*HT *U0*C0 2 2 2 2 2 2 2 2 4 + 2*SIN(AX) *HX *HT *U0 - SIN(AX) *HX *HT *C0 + 2*HX )/(2 4 4 4 2 2 2 2 4 *(SIN(AX) *HT *U0 + 2*SIN(AX) *HX *HT *U0 + HX )) 2 3 2 2 AI22 := (SIN(AX)*HT*( - SIN(AX) *COS(AX)*HT *U0 *C0 2 2 3 2 2 2 - 2*SIN(AX) *HX*HT *U0 + 2*SIN(AX) *HX*HT *U0*C0 2 2 3 + COS(AX)*HX *HT*C0 - 2*HX *U0))/(2 4 4 4 2 2 2 2 4 *(SIN(AX) *HT *U0 + 2*SIN(AX) *HX *HT *U0 + HX )) cleardenot; clear aa,bb,pol; %*********************************************************************** % Example C.4 : Richtmyer, Morton: Difference methods for initial value % problems, &10.3. p.262 coordinates t,x into n,j; unfunc v,w; matrix aa(1,2),bb(2,2); aa(1,1):=(v(n+1)-v(n))/ht-c*(w(j+1/2)-w(j-1/2)+ w(n+1,j+1/2)-w(n+1,j-1/2))/(2*hx)$ aa(1,2):=(w(n+1,j-1/2)-w(n,j-1/2))/ht-c*(v(n+1,j)-v(n+1,j-1)+ v(j)-v(j-1))/(2*hx)$ bb:=ampmat aa; KX*HX AX := ------- 2 H1 := MAT((HX, - SIN(AX)*C*I*HT), 2 2 (C*HT*( - SIN(AX) - 2*SIN(AX)*COS(AX)*I + COS(AX) - 1), 2*HX*( - SIN(AX)*I + COS(AX)))) H0 := MAT((HX,SIN(AX)*C*I*HT), 2 2 (C*HT*(SIN(AX) + 2*SIN(AX)*COS(AX)*I - COS(AX) + 1), 2*HX*( - SIN(AX)*I + COS(AX)))) [ 2 2 2 2 ] [ - SIN(AX) *C *HT + HX 2*SIN(AX)*C*I*HX*HT ] [-------------------------- ----------------------- ] [ 2 2 2 2 2 2 2 2 ] [ SIN(AX) *C *HT + HX SIN(AX) *C *HT + HX ] BB := [ ] [ 2 2 2 2 ] [ 2*SIN(AX)*C*I*HX*HT - SIN(AX) *C *HT + HX ] [ ----------------------- --------------------------] [ 2 2 2 2 2 2 2 2 ] [ SIN(AX) *C *HT + HX SIN(AX) *C *HT + HX ] bb:=denotemat bb; [ AR11 I*AI12] BB := [ ] [I*AI21 AR22 ] pol:=charpol bb; 2 POL := LAM - LAM*(AR11 + AR22) + AI12*AI21 + AR11*AR22 prdenot; 2 2 2 2 - SIN(AX) *C *HT + HX AR11 := -------------------------- 2 2 2 2 SIN(AX) *C *HT + HX 2*SIN(AX)*C*HX*HT AI12 := ----------------------- 2 2 2 2 SIN(AX) *C *HT + HX 2*SIN(AX)*C*HX*HT AI21 := ----------------------- 2 2 2 2 SIN(AX) *C *HT + HX 2 2 2 2 - SIN(AX) *C *HT + HX AR22 := -------------------------- 2 2 2 2 SIN(AX) *C *HT + HX cleardenot; clear aa,bb,pol; %*********************************************************************** % Example C.5: Mazurik: Algoritmy resenia zadaci..., Preprint no.24-85, % AN USSR SO, Inst. teor. i prikl. mechaniky, p.34 coordinates t,x,y into n,m,k; grid uniform,t,x,y; unfunc u1,u2,u3; matrix aa(1,3),bb(3,3); aa(1,1):=(u1(n+1)-u1(n))/ht+c/2*((-u1(m-1)+2*u1(m)-u1(m+1))/hx + (u2(m+1)-u2(m-1))/hx - (u1(k-1)-2*u1(k)+u1(k+1))/hy + (u3(k+1)-u3(k-1))/hy)$ aa(1,2):=(u2(n+1)-u2(n))/ht+c/2*((u1(m+1)-u1(m-1))/hx - (u2(m-1)-2*u2(m)+u2(m+1))/hx)$ aa(1,3):=(u3(n+1)-u3(n))/ht + c/2*((u1(k+1)-u1(k-1))/hy - (u3(k-1)-2*u3(k)+u3(k+1))/hy)$ off prfourmat; bb:=ampmat aa; AX := KX*HX AY := KY*HY BB := MAT(((COS(AX)*C*HT*HY + COS(AY)*C*HX*HT - C*HX*HT - C*HT*HY - SIN(AX)*C*I*HT - SIN(AY)*C*I*HT + HX*HY)/(HX*HY),-------------------,------------------- HX HY ), - SIN(AX)*C*I*HT COS(AX)*C*HT - C*HT + HX (-------------------,--------------------------,0), HX HX - SIN(AY)*C*I*HT COS(AY)*C*HT - C*HT + HY (-------------------,0,--------------------------)) HY HY pol:=charpol bb; 3 2 2 2 POL := (LAM *HX *HY + LAM *HX*HY*( - 2*COS(AX)*C*HT*HY - 2*COS(AY)*C*HX*HT + 2*C*HX*HT + 2*C*HT*HY - 3*HX*HY) + 2 2 2 2 2 2 LAM*(SIN(AY) *C *HX *HT + 3*COS(AX)*COS(AY)*C *HX*HT *HY 2 2 2 2 2 - 3*COS(AX)*C *HX*HT *HY - 2*COS(AX)*C *HT *HY 2 2 2 2 2 + 4*COS(AX)*C*HX*HT*HY + COS(AY) *C *HX *HT 2 2 2 2 2 - 2*COS(AY)*C *HX *HT - 3*COS(AY)*C *HX*HT *HY 2 2 2 2 2 2 + 4*COS(AY)*C*HX *HT*HY + C *HX *HT + 3*C *HX*HT *HY 2 2 2 2 2 + 2*C *HT *HY - 4*C*HX *HT*HY - 4*C*HX*HT*HY 2 2 2 3 3 + 3*HX *HY ) - SIN(AY) *COS(AX)*C *HX*HT 2 3 3 2 2 2 2 + SIN(AY) *C *HX*HT - SIN(AY) *C *HX *HT 2 3 3 3 3 - COS(AX)*COS(AY) *C *HX*HT + 2*COS(AX)*COS(AY)*C *HX*HT 3 3 + 2*COS(AX)*COS(AY)*C *HT *HY 2 2 3 3 - 3*COS(AX)*COS(AY)*C *HX*HT *HY - COS(AX)*C *HX*HT 3 3 2 2 - 2*COS(AX)*C *HT *HY + 3*COS(AX)*C *HX*HT *HY 2 2 2 2 + 2*COS(AX)*C *HT *HY - 2*COS(AX)*C*HX*HT*HY 2 3 3 2 2 2 2 + COS(AY) *C *HX*HT - COS(AY) *C *HX *HT 3 3 3 3 - 2*COS(AY)*C *HX*HT - 2*COS(AY)*C *HT *HY 2 2 2 2 2 + 2*COS(AY)*C *HX *HT + 3*COS(AY)*C *HX*HT *HY 2 3 3 3 3 - 2*COS(AY)*C*HX *HT*HY + C *HX*HT + 2*C *HT *HY 2 2 2 2 2 2 2 2 2 - C *HX *HT - 3*C *HX*HT *HY - 2*C *HT *HY + 2*C*HX *HT*HY 2 2 2 2 2 + 2*C*HX*HT*HY - HX *HY )/(HX *HY ) let cos ax=cos ax2**2-sin ax2**2, cos ay=cos ay2**2-sin ay2**2, sin ax=2*sin ax2*cos ax2, sin ay=2*sin ay2*cos ay2, cos ax2**2=1-sin ax2**2, cos ay2**2=1-sin ay2**2, sin ax2=s1, sin ay2=s2, hx=c*ht/cap1, hy=c*ht/cap2; order s1,s2; pol:=pol; 3 2 2 2 POL := LAM + LAM *(4*S1 *CAP1 + 4*S2 *CAP2 - 3) + LAM*( 2 2 2 2 2 12*S1 *S2 *CAP1*CAP2 + 4*S1 *CAP1 - 8*S1 *CAP1 2 2 2 2 2 2 + 4*S2 *CAP2 - 8*S2 *CAP2 + 3) + 8*S1 *S2 *CAP1 *CAP2 2 2 2 2 2 2 2 + 8*S1 *S2 *CAP1*CAP2 - 12*S1 *S2 *CAP1*CAP2 - 4*S1 *CAP1 2 2 2 2 + 4*S1 *CAP1 - 4*S2 *CAP2 + 4*S2 *CAP2 - 1 clear cos ax,cos ay,sin ax,sin ay,cos ax2**2,cos ay2**2,sin ax2,sin ay2, hx,hy; pol:=complexpol pol; 2 2 If 8*S1 *S2 *CAP1*CAP2*(CAP1 + CAP2) = 0 and 0 = 0 , a root of the polynomial is equal to 1. 3 2 2 2 POL := LAM + LAM *(4*S1 *CAP1 + 4*S2 *CAP2 - 3) + LAM*( 2 2 2 2 2 12*S1 *S2 *CAP1*CAP2 + 4*S1 *CAP1 - 8*S1 *CAP1 2 2 2 2 2 2 + 4*S2 *CAP2 - 8*S2 *CAP2 + 3) + 8*S1 *S2 *CAP1 *CAP2 2 2 2 2 2 2 2 + 8*S1 *S2 *CAP1*CAP2 - 12*S1 *S2 *CAP1*CAP2 - 4*S1 *CAP1 2 2 2 2 + 4*S1 *CAP1 - 4*S2 *CAP2 + 4*S2 *CAP2 - 1 pol1:=hurw pol; 3 2 2 2 POL1 := 8*LAM *S1 *S2 *CAP1*CAP2*(CAP1 + CAP2) + 8*LAM *( 2 2 2 2 2 2 - 3*S1 *S2 *CAP1 *CAP2 - 3*S1 *S2 *CAP1*CAP2 2 2 2 2 2 2 + 3*S1 *S2 *CAP1*CAP2 + S1 *CAP1 + S2 *CAP2 ) + 8*LAM*( 2 2 2 2 2 2 3*S1 *S2 *CAP1 *CAP2 + 3*S1 *S2 *CAP1*CAP2 2 2 2 2 2 - 6*S1 *S2 *CAP1*CAP2 - 2*S1 *CAP1 + 2*S1 *CAP1 2 2 2 2 2 2 - 2*S2 *CAP2 + 2*S2 *CAP2) + 8*( - S1 *S2 *CAP1 *CAP2 2 2 2 2 2 2 2 - S1 *S2 *CAP1*CAP2 + 3*S1 *S2 *CAP1*CAP2 + S1 *CAP1 2 2 2 2 - 2*S1 *CAP1 + S2 *CAP2 - 2*S2 *CAP2 + 1) denotid cp; pol:=denotepol pol; 3 2 POL := LAM + LAM *CPR02 + LAM*CPR01 + CPR00 prdenot; 2 2 2 2 2 2 CPR00 := 8*S1 *S2 *CAP1 *CAP2 + 8*S1 *S2 *CAP1*CAP2 2 2 2 2 2 - 12*S1 *S2 *CAP1*CAP2 - 4*S1 *CAP1 + 4*S1 *CAP1 2 2 2 - 4*S2 *CAP2 + 4*S2 *CAP2 - 1 2 2 2 2 2 CPR01 := 12*S1 *S2 *CAP1*CAP2 + 4*S1 *CAP1 - 8*S1 *CAP1 2 2 2 + 4*S2 *CAP2 - 8*S2 *CAP2 + 3 2 2 CPR02 := 4*S1 *CAP1 + 4*S2 *CAP2 - 3 cleardenot; clear aa,bb,pol,pol1; %*********************************************************************** % Example C.6 : Lax-Wendrov (V. Ganzha) coordinates t,x,y into n,m,k; grid uniform,t,x,y; let cos ax**2=1-sin ax**2, cos ay**2=1-sin ay**2; unfunc u1,u2,u3,u4; matrix aa(1,4),bb(4,4); aa(1,1):=4*(u1(n+1)-u1(n))/ht+ (w*(u1(m+2)-u1(m-2)+u1(m+1,k+1)+u1(m+1,k-1)- u1(m-1,k+1)-u1(m-1,k-1))+p*(u2(m+2)-u2(m-2)+u2(m+1,k+1)+ u2(m+1,k-1)-u2(m-1,k+1)-u2(m-1,k-1))+ v*(u1(m+1,k+1)+u1(m-1,k+1)- u1(m+1,k-1)-u1(m-1,k-1)+u1(k+2)-u1(k-2))+p*(u3(m+1,k+1)+ u3(m-1,k+1)-u3(m+1,k-1)-u3(m-1,k-1)+u3(k+2)-u3(k-2)))/hx+ht* (2*w**2*(-u1(m+2)+2*u1(m)-u1(m-2))+4*w*p*(-u2(m+2)+2*u2(m)- u2(m-2))+2*(-u4(m+2)+2*u4(m)-u4(m-2))+2*v**2*(-u1(k+2)+ 2*u1(k)-u1(k-2))+4*v*p*(u3(k+2)+2*u3(k)-u3(k-2))+2*(-u4(k+2)+ 2*u4(k)-u4(k-2))+4*w*v*(-u1(m+1,k+1)+u1(m+1,k-1)+u1(m-1,k+1)- u1(m-1,k-1))+4*p*v*(-u2(m+1,k+1)+u2(m+1,k-1)+u2(m-1,k+1)- u2(m-1,k-1))+4*w*p*(-u3(m+1,k+1)+u3(m+1,k-1)+u3(m-1,k+1)- u3(m-1,k-1)))/hx/hx$ aa(1,2):=4*p*(u2(n+1)-u2(n))/ht+ (w*p*(u2(m+2)-u2(m-2)+u2(m+1,k+1)+ u2(m+1,k-1)-u2(m-1,k+1)-u2(m-1,k-1))+u4(m+2)-u4(m-2)+ u4(m+1,k+1)+ u4(m+1,k-1)-u4(m-1,k+1)-u4(m-1,k-1)+ p*v*(u2(m+1,k+1)+u2(m-1,k+1)+ u2(k+2)-u2(k-2)-u2(m+1,k-1)-u2(m-1,k-1)))/hx+ht*(2*w**2*p* (-u2(m+2)+2*u2(m)-u2(m-2))+2*p*c**2*(-u2(m+2)+2*u2(m)-u2(m-2)) +4*w*(-u4(m+2)+2*u4(m)-u4(m-2))+2*p*v**2*(-u2(k+2)+2*u2(k)- u2(k-2))+4*w*p*v*(-u2(m+1,k+1)+u2(m+1,k-1)+u2(m-1,k+1)- u2(m-1,k-1))+2*p*c**2*(-u3(m+1,k+1)+u3(m+1,k-1)+u3(m-1,k+1) -u3(m-1,k-1))+4*v*(-u4(m+1,k+1)+u4(m+1,k-1)+u4(m-1,k+1)- u4(m-1,k-1)))/hx/hx$ aa(1,3):=4*p*(u3(n+1)-u3(n))/ht+(w*p*(u3(m+2)-u3(m-2)+u3(m+1,k+1)+ u3(m+1,k-1)-u3(m-1,k+1)-u3(m-1,k-1))+u4(k+2)-u4(k-2)+ u4(m+1,k+1)-u4(m+1,k-1)+u4(m-1,k+1)-u4(m-1,k-1)+ p*v*(u3(m+1,k+1)+u3(m-1,k+1)+u3(k+2)-u3(k-2)-u3(m+1,k-1)- u3(m-1,k-1)))/hx+ht*(2*w**2*p*(-u3(m+2)+2*u3(m)-u3(m-2))+ 2*p*c**2*(-u3(k+2)+2*u3(k)-u3(k-2))+4*v*(-u4(k+2)+ 2*u4(k)-u4(k-2))+2*p*v**2*(-u3(k+2)+2*u3(k)-u3(k-2))+ 4*w*p*v*(-u3(m+1,k+1)+u3(m+1,k-1)+u3(m-1,k+1)- u3(m-1,k-1))+2*p*c**2*(-u2(m+1,k+1)+u2(m+1,k-1)+ u2(m-1,k+1)-u2(m-1,k-1))+4*w*(u4(m+1,k+1)+u4(m+1,k-1)+ u4(m-1,k+1)-u4(m-1,k-1)))/hx/hx$ aa(1,4):=4*(u4(n+1)-u4(n))/ht+(p*c**2*(u2(m+2)-u2(m-2)+u2(m+1,k+1)+ u2(m+1,k-1)-u2(m-1,k+1)-u2(m-1,k-1))+w*(u4(m+2)- u4(m-2)+u4(m+1,k+1)+u4(m+1,k-1)-u4(m-1,k+1)-u4(m-1,k-1))+ +p*c**2*(u3(m+1,k+1)+u3(m-1,k+1)-u3(m+1,k-1)- u3(m-1,k-1)+u3(k+2)-u3(k-2))+v*(u4(m+1,k+1)+u4(m-1,k+1)- u4(m+1,k-1)-u4(m-1,k-1)+u4(k+2)-u4(k-2)))/hx+ht* (2*w**2*(-u4(m+2)+2*u4(m)-u4(m-2))+4*w*p*c**2*(-u2(m+2)+ 2*u2(m)-u2(m-2))+2*c**2*(-u4(m+2)+2*u4(m)-u4(m-2))+ 4*p*v*c**2*(-u3(k+2)+2*u3(k)-u3(k-2))+2*c**2*(-u4(k+2)+ 2*u4(k)-u4(k-2))+2*v**2*(-u4(k+2)+2*u4(k)-u4(k-2))+ 4*p*v*c**2*(-u2(m+1,k+1)+u2(m+1,k-1)+u2(m-1,k+1)- u2(m-1,k-1))+4*w*p*c**2*(-u3(m+1,k+1)+u3(m+1,k-1)+ u3(m-1,k+1)-u3(m-1,k-1))+4*w*v*(-u4(m+1,k+1)+ u4(m+1,k-1)+u4(m-1,k+1)-u4(m-1,k-1)))/hx/hx$ bb:=ampmat aa; AX := KX*HX AY := KY*HY 2 2 2 2 BB := MAT((( - 2*SIN(AX) *W *HT - 4*SIN(AX)*SIN(AY)*V*W*HT - SIN(AX)*COS(AX)*I*W*HX*HT - SIN(AX)*COS(AY)*I*W*HX*HT 2 2 2 - 2*SIN(AY) *V *HT - SIN(AY)*COS(AX)*I*V*HX*HT 2 2 - SIN(AY)*COS(AY)*I*V*HX*HT + HX )/HX ,(SIN(AX)*P*HT*( - 4*SIN(AX)*W*HT - 4*SIN(AY)*V*HT - COS(AX)*I*HX 2 - COS(AY)*I*HX))/HX ,(P*HT*( - 4*SIN(AX)*SIN(AY)*W*HT - SIN(AY)*COS(AX)*I*HX - 4*SIN(AY)*COS(AY)*I*V*HT - SIN(AY)*COS(AY)*I*HX 2 2 2 2 - (2*HT )*(SIN(AX) + SIN(AY) ) - 2*V*HT))/HX ,----------------------------------), 2 HX 2 2 2 2 2 2 (0,( - 2*SIN(AX) *C *HT - 2*SIN(AX) *W *HT 2 - 4*SIN(AX)*SIN(AY)*V*W*HT - SIN(AX)*COS(AX)*I*W*HX*HT 2 2 2 - SIN(AX)*COS(AY)*I*W*HX*HT - 2*SIN(AY) *V *HT - SIN(AY)*COS(AX)*I*V*HX*HT 2 2 - SIN(AY)*COS(AY)*I*V*HX*HT + HX )/HX , 2 2 - 2*SIN(AX)*SIN(AY)*C *HT -----------------------------,(SIN(AX)*HT*( 2 HX - 4*SIN(AX)*W*HT - 4*SIN(AY)*V*HT - COS(AX)*I*HX 2 - COS(AY)*I*HX))/(P*HX )), 2 2 - 2*SIN(AX)*SIN(AY)*C *HT 2 2 2 (0,-----------------------------,( - 2*SIN(AX) *W *HT 2 HX 2 - 4*SIN(AX)*SIN(AY)*V*W*HT - SIN(AX)*COS(AX)*I*W*HX*HT 2 2 2 - SIN(AX)*COS(AY)*I*W*HX*HT - 2*SIN(AY) *C *HT 2 2 2 - 2*SIN(AY) *V *HT - SIN(AY)*COS(AX)*I*V*HX*HT 2 2 - SIN(AY)*COS(AY)*I*V*HX*HT + HX )/HX ,(HT*( - 2*SIN(AX)*SIN(AY)*W*HT - 2*SIN(AX)*COS(AY)*I*W*HT 2 - 4*SIN(AY) *V*HT - 2*SIN(AY)*COS(AX)*I*W*HT - SIN(AY)*COS(AX)*I*HX - SIN(AY)*COS(AY)*I*HX 2 - 2*COS(AX)*COS(AY)*W*HT))/(P*HX )), 2 (0,(SIN(AX)*C *P*HT*( - 4*SIN(AX)*W*HT - 4*SIN(AY)*V*HT 2 2 - COS(AX)*I*HX - COS(AY)*I*HX))/HX ,(SIN(AY)*C *P *HT*( - 4*SIN(AX)*W*HT - 4*SIN(AY)*V*HT - COS(AX)*I*HX 2 2 2 2 - COS(AY)*I*HX))/HX ,( - 2*SIN(AX) *C *HT 2 2 2 2 - 2*SIN(AX) *W *HT - 4*SIN(AX)*SIN(AY)*V*W*HT - SIN(AX)*COS(AX)*I*W*HX*HT 2 2 2 - SIN(AX)*COS(AY)*I*W*HX*HT - 2*SIN(AY) *C *HT 2 2 2 - 2*SIN(AY) *V *HT - SIN(AY)*COS(AX)*I*V*HX*HT 2 2 - SIN(AY)*COS(AY)*I*V*HX*HT + HX )/HX )) let sin(ax)=s1, cos(ax)=c1, sin(ay)=s2, cos(ay)=c2, w=k1*hx/ht, v=k2*hx/ht, c=k3*hx/ht, ht=r1*hx; denotid a; bb:=denotemat bb; [I*AI11 + AR11 I*AI12 + AR12 I*AI13 + AR13 AR14 ] [ ] [ 0 I*AI22 + AR22 AR23 I*AI24 + AR24] BB := [ ] [ 0 AR32 I*AI33 + AR33 I*AI34 + AR34] [ ] [ 0 I*AI42 + AR42 I*AI43 + AR43 I*AI44 + AR44] clear sin ax,cos ax,sin ay,cos ay,w,v,c,ht; pol:=charpol bb; 4 3 POL := LAM - LAM *(I*AI11 + I*AI22 + I*AI33 + I*AI44 + AR11 + AR22 2 + AR33 + AR44) + LAM *(I*AR11*AI22 + I*AR11*AI33 + I*AR11*AI44 + I*AI11*AR22 + I*AI11*AR33 + I*AI11*AR44 + I*AR22*AI33 + I*AR22*AI44 + I*AI22*AR33 + I*AI22*AR44 - I*AR24*AI42 - I*AI24*AR42 + I*AR33*AI44 + I*AI33*AR44 - I*AR34*AI43 - I*AI34*AR43 + AR11*AR22 + AR11*AR33 + AR11*AR44 - AI11*AI22 - AI11*AI33 - AI11*AI44 + AR22*AR33 + AR22*AR44 - AI22*AI33 - AI22*AI44 - AR23*AR32 - AR24*AR42 + AI24*AI42 + AR33*AR44 - AI33*AI44 - AR34*AR43 + AI34*AI43) + LAM*( - I*AR11*AR22*AI33 - I*AR11*AR22*AI44 - I*AR11*AI22*AR33 - I*AR11*AI22*AR44 + I*AR11*AR24*AI42 + I*AR11*AI24*AR42 - I*AR11*AR33*AI44 - I*AR11*AI33*AR44 + I*AR11*AR34*AI43 + I*AR11*AI34*AR43 - I*AI11*AR22*AR33 - I*AI11*AR22*AR44 + I*AI11*AI22*AI33 + I*AI11*AI22*AI44 + I*AI11*AR23*AR32 + I*AI11*AR24*AR42 - I*AI11*AI24*AI42 - I*AI11*AR33*AR44 + I*AI11*AI33*AI44 + I*AI11*AR34*AR43 - I*AI11*AI34*AI43 - I*AR22*AR33*AI44 - I*AR22*AI33*AR44 + I*AR22*AR34*AI43 + I*AR22*AI34*AR43 - I*AI22*AR33*AR44 + I*AI22*AI33*AI44 + I*AI22*AR34*AR43 - I*AI22*AI34*AI43 + I*AR23*AR32*AI44 - I*AR23*AR34*AI42 - I*AR23*AI34*AR42 - I*AR24*AR32*AI43 + I*AR24*AR33*AI42 + I*AR24*AI33*AR42 - I*AI24*AR32*AR43 + I*AI24*AR33*AR42 - I*AI24*AI33*AI42 - AR11*AR22*AR33 - AR11*AR22*AR44 + AR11*AI22*AI33 + AR11*AI22*AI44 + AR11*AR23*AR32 + AR11*AR24*AR42 - AR11*AI24*AI42 - AR11*AR33*AR44 + AR11*AI33*AI44 + AR11*AR34*AR43 - AR11*AI34*AI43 + AI11*AR22*AI33 + AI11*AR22*AI44 + AI11*AI22*AR33 + AI11*AI22*AR44 - AI11*AR24*AI42 - AI11*AI24*AR42 + AI11*AR33*AI44 + AI11*AI33*AR44 - AI11*AR34*AI43 - AI11*AI34*AR43 - AR22*AR33*AR44 + AR22*AI33*AI44 + AR22*AR34*AR43 - AR22*AI34*AI43 + AI22*AR33*AI44 + AI22*AI33*AR44 - AI22*AR34*AI43 - AI22*AI34*AR43 + AR23*AR32*AR44 - AR23*AR34*AR42 + AR23*AI34*AI42 - AR24*AR32*AR43 + AR24*AR33*AR42 - AR24*AI33*AI42 + AI24*AR32*AI43 - AI24*AR33*AI42 - AI24*AI33*AR42) + I*AR11*AR22*AR33*AI44 + I*AR11*AR22*AI33*AR44 - I*AR11*AR22*AR34*AI43 - I*AR11*AR22*AI34*AR43 + I*AR11*AI22*AR33*AR44 - I*AR11*AI22*AI33*AI44 - I*AR11*AI22*AR34*AR43 + I*AR11*AI22*AI34*AI43 - I*AR11*AR23*AR32*AI44 + I*AR11*AR23*AR34*AI42 + I*AR11*AR23*AI34*AR42 + I*AR11*AR24*AR32*AI43 - I*AR11*AR24*AR33*AI42 - I*AR11*AR24*AI33*AR42 + I*AR11*AI24*AR32*AR43 - I*AR11*AI24*AR33*AR42 + I*AR11*AI24*AI33*AI42 + I*AI11*AR22*AR33*AR44 - I*AI11*AR22*AI33*AI44 - I*AI11*AR22*AR34*AR43 + I*AI11*AR22*AI34*AI43 - I*AI11*AI22*AR33*AI44 - I*AI11*AI22*AI33*AR44 + I*AI11*AI22*AR34*AI43 + I*AI11*AI22*AI34*AR43 - I*AI11*AR23*AR32*AR44 + I*AI11*AR23*AR34*AR42 - I*AI11*AR23*AI34*AI42 + I*AI11*AR24*AR32*AR43 - I*AI11*AR24*AR33*AR42 + I*AI11*AR24*AI33*AI42 - I*AI11*AI24*AR32*AI43 + I*AI11*AI24*AR33*AI42 + I*AI11*AI24*AI33*AR42 + AR11*AR22*AR33*AR44 - AR11*AR22*AI33*AI44 - AR11*AR22*AR34*AR43 + AR11*AR22*AI34*AI43 - AR11*AI22*AR33*AI44 - AR11*AI22*AI33*AR44 + AR11*AI22*AR34*AI43 + AR11*AI22*AI34*AR43 - AR11*AR23*AR32*AR44 + AR11*AR23*AR34*AR42 - AR11*AR23*AI34*AI42 + AR11*AR24*AR32*AR43 - AR11*AR24*AR33*AR42 + AR11*AR24*AI33*AI42 - AR11*AI24*AR32*AI43 + AR11*AI24*AR33*AI42 + AR11*AI24*AI33*AR42 - AI11*AR22*AR33*AI44 - AI11*AR22*AI33*AR44 + AI11*AR22*AR34*AI43 + AI11*AR22*AI34*AR43 - AI11*AI22*AR33*AR44 + AI11*AI22*AI33*AI44 + AI11*AI22*AR34*AR43 - AI11*AI22*AI34*AI43 + AI11*AR23*AR32*AI44 - AI11*AR23*AR34*AI42 - AI11*AR23*AI34*AR42 - AI11*AR24*AR32*AI43 + AI11*AR24*AR33*AI42 + AI11*AR24*AI33*AR42 - AI11*AI24*AR32*AR43 + AI11*AI24*AR33*AR42 - AI11*AI24*AI33*AI42 denotid cp; pol:=denotepol pol; 4 3 2 POL := LAM + LAM *(I*CPI03 + CPR03) + LAM *(I*CPI02 + CPR02) + LAM*(I*CPI01 + CPR01) + I*CPI00 + CPR00 pol:=complexpol pol; If CPR00 + CPR01 + CPR02 + CPR03 + 1 = 0 and CPI00 + CPI01 + CPI02 + CPI03 = 0 , a root of the polynomial is equal to 1. 8 7 6 2 2 POL := LAM + 2*LAM *CPR03 + LAM *(2*CPR02 + CPR03 + CPI03 ) 5 4 + 2*LAM *(CPR01 + CPR02*CPR03 + CPI02*CPI03) + LAM 2 2 *(2*CPR00 + 2*CPR01*CPR03 + 2*CPI01*CPI03 + CPR02 + CPI02 ) + 3 2*LAM *(CPR00*CPR03 + CPI00*CPI03 + CPR01*CPR02 + CPI01*CPI02) 2 2 2 + LAM *(2*CPR00*CPR02 + 2*CPI00*CPI02 + CPR01 + CPI01 ) 2 2 + 2*LAM*(CPR00*CPR01 + CPI00*CPI01) + CPR00 + CPI00 denotid rp; pol:=denotepol pol; 8 7 6 5 4 POL := LAM + LAM *RPR07 + LAM *RPR06 + LAM *RPR05 + LAM *RPR04 3 2 + LAM *RPR03 + LAM *RPR02 + LAM*RPR01 + RPR00 prdenot; 2 2 2 2 AR11 := - 2*S1 *K1 - 4*S1*S2*K1*K2 - 2*S2 *K2 + 1 AI11 := - (S1*C1*K1 + S1*C2*K1 + S2*C1*K2 + S2*C2*K2) AR12 := - (4*S1*P*R1)*(S1*K1 + S2*K2) AI12 := - (S1*P*R1)*(C1 + C2) AR13 := 2*P*R1*( - 2*S1*S2*K1 - K2) AI13 := S2*P*R1*( - C1 - 4*C2*K2 - C2) 2 2 2 AR14 := - (2*R1 )*(S1 + S2 ) 2 2 2 2 2 2 AR22 := - 2*S1 *K1 - 2*S1 *K3 - 4*S1*S2*K1*K2 - 2*S2 *K2 + 1 AI22 := - (S1*C1*K1 + S1*C2*K1 + S2*C1*K2 + S2*C2*K2) 2 AR23 := - 2*S1*S2*K3 - (4*S1*R1)*(S1*K1 + S2*K2) AR24 := ------------------------------ P - (S1*R1)*(C1 + C2) AI24 := ---------------------- P 2 AR32 := - 2*S1*S2*K3 2 2 2 2 2 2 AR33 := - 2*S1 *K1 - 4*S1*S2*K1*K2 - 2*S2 *K2 - 2*S2 *K3 + 1 AI33 := - (S1*C1*K1 + S1*C2*K1 + S2*C1*K2 + S2*C2*K2) 2 2*R1*( - S1*S2*K1 - 2*S2 *K2 - C1*C2*K1) AR34 := ------------------------------------------ P R1*( - 2*S1*C2*K1 - 2*S2*C1*K1 - S2*C1 - S2*C2) AI34 := ------------------------------------------------- P 2 - (4*S1*P*K3 )*(S1*K1 + S2*K2) AR42 := --------------------------------- R1 2 - (S1*P*K3 )*(C1 + C2) AI42 := ------------------------- R1 2 - (4*S2*P*K3 )*(S1*K1 + S2*K2) AR43 := --------------------------------- R1 2 - (S2*P*K3 )*(C1 + C2) AI43 := ------------------------- R1 AR44 := 2 2 2 2 2 2 2 2 - 2*S1 *K1 - 2*S1 *K3 - 4*S1*S2*K1*K2 - 2*S2 *K2 - 2*S2 *K3 + 1 AI44 := - (S1*C1*K1 + S1*C2*K1 + S2*C1*K2 + S2*C2*K2) CPR00 := AR11*AR22*AR33*AR44 - AR11*AR22*AI33*AI44 - AR11*AR22*AR34*AR43 + AR11*AR22*AI34*AI43 - AR11*AI22*AR33*AI44 - AR11*AI22*AI33*AR44 + AR11*AI22*AR34*AI43 + AR11*AI22*AI34*AR43 - AR11*AR23*AR32*AR44 + AR11*AR23*AR34*AR42 - AR11*AR23*AI34*AI42 + AR11*AR24*AR32*AR43 - AR11*AR24*AR33*AR42 + AR11*AR24*AI33*AI42 - AR11*AI24*AR32*AI43 + AR11*AI24*AR33*AI42 + AR11*AI24*AI33*AR42 - AI11*AR22*AR33*AI44 - AI11*AR22*AI33*AR44 + AI11*AR22*AR34*AI43 + AI11*AR22*AI34*AR43 - AI11*AI22*AR33*AR44 + AI11*AI22*AI33*AI44 + AI11*AI22*AR34*AR43 - AI11*AI22*AI34*AI43 + AI11*AR23*AR32*AI44 - AI11*AR23*AR34*AI42 - AI11*AR23*AI34*AR42 - AI11*AR24*AR32*AI43 + AI11*AR24*AR33*AI42 + AI11*AR24*AI33*AR42 - AI11*AI24*AR32*AR43 + AI11*AI24*AR33*AR42 - AI11*AI24*AI33*AI42 CPI00 := AR11*AR22*AR33*AI44 + AR11*AR22*AI33*AR44 - AR11*AR22*AR34*AI43 - AR11*AR22*AI34*AR43 + AR11*AI22*AR33*AR44 - AR11*AI22*AI33*AI44 - AR11*AI22*AR34*AR43 + AR11*AI22*AI34*AI43 - AR11*AR23*AR32*AI44 + AR11*AR23*AR34*AI42 + AR11*AR23*AI34*AR42 + AR11*AR24*AR32*AI43 - AR11*AR24*AR33*AI42 - AR11*AR24*AI33*AR42 + AR11*AI24*AR32*AR43 - AR11*AI24*AR33*AR42 + AR11*AI24*AI33*AI42 + AI11*AR22*AR33*AR44 - AI11*AR22*AI33*AI44 - AI11*AR22*AR34*AR43 + AI11*AR22*AI34*AI43 - AI11*AI22*AR33*AI44 - AI11*AI22*AI33*AR44 + AI11*AI22*AR34*AI43 + AI11*AI22*AI34*AR43 - AI11*AR23*AR32*AR44 + AI11*AR23*AR34*AR42 - AI11*AR23*AI34*AI42 + AI11*AR24*AR32*AR43 - AI11*AR24*AR33*AR42 + AI11*AR24*AI33*AI42 - AI11*AI24*AR32*AI43 + AI11*AI24*AR33*AI42 + AI11*AI24*AI33*AR42 CPR01 := - AR11*AR22*AR33 - AR11*AR22*AR44 + AR11*AI22*AI33 + AR11*AI22*AI44 + AR11*AR23*AR32 + AR11*AR24*AR42 - AR11*AI24*AI42 - AR11*AR33*AR44 + AR11*AI33*AI44 + AR11*AR34*AR43 - AR11*AI34*AI43 + AI11*AR22*AI33 + AI11*AR22*AI44 + AI11*AI22*AR33 + AI11*AI22*AR44 - AI11*AR24*AI42 - AI11*AI24*AR42 + AI11*AR33*AI44 + AI11*AI33*AR44 - AI11*AR34*AI43 - AI11*AI34*AR43 - AR22*AR33*AR44 + AR22*AI33*AI44 + AR22*AR34*AR43 - AR22*AI34*AI43 + AI22*AR33*AI44 + AI22*AI33*AR44 - AI22*AR34*AI43 - AI22*AI34*AR43 + AR23*AR32*AR44 - AR23*AR34*AR42 + AR23*AI34*AI42 - AR24*AR32*AR43 + AR24*AR33*AR42 - AR24*AI33*AI42 + AI24*AR32*AI43 - AI24*AR33*AI42 - AI24*AI33*AR42 CPI01 := - AR11*AR22*AI33 - AR11*AR22*AI44 - AR11*AI22*AR33 - AR11*AI22*AR44 + AR11*AR24*AI42 + AR11*AI24*AR42 - AR11*AR33*AI44 - AR11*AI33*AR44 + AR11*AR34*AI43 + AR11*AI34*AR43 - AI11*AR22*AR33 - AI11*AR22*AR44 + AI11*AI22*AI33 + AI11*AI22*AI44 + AI11*AR23*AR32 + AI11*AR24*AR42 - AI11*AI24*AI42 - AI11*AR33*AR44 + AI11*AI33*AI44 + AI11*AR34*AR43 - AI11*AI34*AI43 - AR22*AR33*AI44 - AR22*AI33*AR44 + AR22*AR34*AI43 + AR22*AI34*AR43 - AI22*AR33*AR44 + AI22*AI33*AI44 + AI22*AR34*AR43 - AI22*AI34*AI43 + AR23*AR32*AI44 - AR23*AR34*AI42 - AR23*AI34*AR42 - AR24*AR32*AI43 + AR24*AR33*AI42 + AR24*AI33*AR42 - AI24*AR32*AR43 + AI24*AR33*AR42 - AI24*AI33*AI42 CPR02 := AR11*AR22 + AR11*AR33 + AR11*AR44 - AI11*AI22 - AI11*AI33 - AI11*AI44 + AR22*AR33 + AR22*AR44 - AI22*AI33 - AI22*AI44 - AR23*AR32 - AR24*AR42 + AI24*AI42 + AR33*AR44 - AI33*AI44 - AR34*AR43 + AI34*AI43 CPI02 := AR11*AI22 + AR11*AI33 + AR11*AI44 + AI11*AR22 + AI11*AR33 + AI11*AR44 + AR22*AI33 + AR22*AI44 + AI22*AR33 + AI22*AR44 - AR24*AI42 - AI24*AR42 + AR33*AI44 + AI33*AR44 - AR34*AI43 - AI34*AR43 CPR03 := - (AR11 + AR22 + AR33 + AR44) CPI03 := - (AI11 + AI22 + AI33 + AI44) 2 2 RPR00 := CPR00 + CPI00 RPR01 := 2*(CPR00*CPR01 + CPI00*CPI01) 2 2 RPR02 := 2*CPR00*CPR02 + 2*CPI00*CPI02 + CPR01 + CPI01 RPR03 := 2*(CPR00*CPR03 + CPI00*CPI03 + CPR01*CPR02 + CPI01*CPI02) 2 2 RPR04 := 2*CPR00 + 2*CPR01*CPR03 + 2*CPI01*CPI03 + CPR02 + CPI02 RPR05 := 2*(CPR01 + CPR02*CPR03 + CPI02*CPI03) 2 2 RPR06 := 2*CPR02 + CPR03 + CPI03 RPR07 := 2*CPR03 cleardenot; clear aa,bb,pol; %*********************************************************************** %***** ***** %***** T e s t Examples --- Module H U R W P ***** %***** ***** %*********************************************************************** % Example H.1 x0:=lam-1; X0 := LAM - 1 x1:=lam-(ar+i*ai); X1 := LAM - (I*AI + AR) x2:=lam-(br+i*bi); X2 := LAM - (I*BI + BR) x3:=lam-(cr+i*ci); X3 := LAM - (I*CI + CR) hurwitzp x1; Necessary and sufficient conditions are: - AR > 0 COND % Example H.2 x:=hurw(x0*x1); X := 2*LAM*( - I*AI - AR + 1) + 2*(I*AI + AR + 1) hurwitzp x; Necessary and sufficient conditions are: 2 2 4*( - AR - AI + 1) > 0 COND % Example H.3 x:=(x1*x2); 2 X := LAM - LAM*(I*AI + I*BI + AR + BR) + I*AR*BI + I*AI*BR + AR*BR - AI*BI hurwitzp x; Necessary and sufficient conditions are: - (AR + BR) > 0 2 2 2 2 AR*BR*(AR + 2*AR*BR + AI - 2*AI*BI + BR + BI ) > 0 COND % Example H.4 x:=(x1*x2*x3); 3 2 X := LAM - LAM *(I*AI + I*BI + I*CI + AR + BR + CR) + LAM*(I*AR*BI + I*AR*CI + I*AI*BR + I*AI*CR + I*BR*CI + I*BI*CR + AR*BR + AR*CR - AI*BI - AI*CI + BR*CR - BI*CI) - I*AR*BR*CI - I*AR*BI*CR - I*AI*BR*CR + I*AI*BI*CI - AR*BR*CR + AR*BI*CI + AI*BR*CI + AI*BI*CR hurwitzp x; Necessary and sufficient conditions are: - (AR + BR + CR) > 0 3 3 2 2 2 2 2 2 AR *BR + AR *CR + 2*AR *BR + 4*AR *BR*CR + 2*AR *CR + AR*AI *BR 2 3 2 + AR*AI *CR - 2*AR*AI*BR*BI - 2*AR*AI*CR*CI + AR*BR + 4*AR*BR *CR 2 2 3 2 3 2 2 + AR*BR*BI + 4*AR*BR*CR + AR*CR + AR*CR*CI + BR *CR + 2*BR *CR 2 3 2 + BR*BI *CR - 2*BR*BI*CR*CI + BR*CR + BR*CR*CI > 0 4 2 4 4 2 4 4 2 AR*BR*CR*( - AR *BR - 2*AR *BR*CR - AR *BI + 2*AR *BI*CI - AR *CR 4 2 3 3 3 2 3 2 - AR *CI - 2*AR *BR - 6*AR *BR *CR - 2*AR *BR*BI 3 3 2 3 2 + 4*AR *BR*BI*CI - 6*AR *BR*CR - 2*AR *BR*CI 3 2 3 3 3 3 2 - 2*AR *BI *CR + 4*AR *BI*CR*CI - 2*AR *CR - 2*AR *CR*CI 2 2 2 2 2 2 2 2 - 2*AR *AI *BR - 4*AR *AI *BR*CR - 2*AR *AI *BI 2 2 2 2 2 2 2 2 + 4*AR *AI *BI*CI - 2*AR *AI *CR - 2*AR *AI *CI 2 2 2 2 2 + 2*AR *AI*BR *BI + 2*AR *AI*BR *CI + 4*AR *AI*BR*BI*CR 2 2 3 2 2 + 4*AR *AI*BR*CR*CI + 2*AR *AI*BI - 2*AR *AI*BI *CI 2 2 2 2 2 2 + 2*AR *AI*BI*CR - 2*AR *AI*BI*CI + 2*AR *AI*CR *CI 2 3 2 4 2 3 2 2 2 + 2*AR *AI*CI - AR *BR - 6*AR *BR *CR - 2*AR *BR *BI 2 2 2 2 2 2 2 2 + 2*AR *BR *BI*CI - 10*AR *BR *CR - 2*AR *BR *CI 2 2 2 2 3 - 6*AR *BR*BI *CR + 8*AR *BR*BI*CR*CI - 6*AR *BR*CR 2 2 2 4 2 3 2 2 2 - 6*AR *BR*CR*CI - AR *BI + 2*AR *BI *CI - 2*AR *BI *CR 2 2 2 2 2 2 3 2 4 - 2*AR *BI *CI + 2*AR *BI*CR *CI + 2*AR *BI*CI - AR *CR 2 2 2 2 4 2 3 2 2 - 2*AR *CR *CI - AR *CI - 2*AR*AI *BR - 6*AR*AI *BR *CR 2 2 2 2 2 - 2*AR*AI *BR*BI + 4*AR*AI *BR*BI*CI - 6*AR*AI *BR*CR 2 2 2 2 2 - 2*AR*AI *BR*CI - 2*AR*AI *BI *CR + 4*AR*AI *BI*CR*CI 2 3 2 2 3 - 2*AR*AI *CR - 2*AR*AI *CR*CI + 4*AR*AI*BR *CI 2 2 + 4*AR*AI*BR *BI*CR + 8*AR*AI*BR *CR*CI 2 2 + 4*AR*AI*BR*BI *CI + 8*AR*AI*BR*BI*CR 2 2 3 - 8*AR*AI*BR*BI*CI + 4*AR*AI*BR*CR *CI + 4*AR*AI*BR*CI 3 2 3 + 4*AR*AI*BI *CR - 8*AR*AI*BI *CR*CI + 4*AR*AI*BI*CR 2 4 3 2 + 4*AR*AI*BI*CR*CI - 2*AR*BR *CR - 6*AR*BR *CR 3 2 2 2 2 - 2*AR*BR *CI - 4*AR*BR *BI *CR + 4*AR*BR *BI*CR*CI 2 3 2 2 2 2 - 6*AR*BR *CR - 6*AR*BR *CR*CI - 6*AR*BR*BI *CR 2 2 2 3 - 2*AR*BR*BI *CI + 4*AR*BR*BI*CR *CI + 4*AR*BR*BI*CI 4 2 2 4 - 2*AR*BR*CR - 4*AR*BR*CR *CI - 2*AR*BR*CI 4 3 2 3 - 2*AR*BI *CR + 4*AR*BI *CR*CI - 2*AR*BI *CR 2 2 4 2 4 4 2 - 2*AR*BI *CR*CI - AI *BR - 2*AI *BR*CR - AI *BI 4 4 2 4 2 3 2 + 2*AI *BI*CI - AI *CR - AI *CI + 2*AI *BR *BI 3 2 3 3 + 2*AI *BR *CI + 4*AI *BR*BI*CR + 4*AI *BR*CR*CI 3 3 3 2 3 2 3 2 + 2*AI *BI - 2*AI *BI *CI + 2*AI *BI*CR - 2*AI *BI*CI 3 2 3 3 2 4 2 3 + 2*AI *CR *CI + 2*AI *CI - AI *BR - 2*AI *BR *CR 2 2 2 2 2 2 2 2 - 2*AI *BR *BI - 2*AI *BR *BI*CI - 2*AI *BR *CR 2 2 2 2 2 2 - 2*AI *BR *CI - 2*AI *BR*BI *CR - 8*AI *BR*BI*CR*CI 2 3 2 2 2 4 2 3 - 2*AI *BR*CR - 2*AI *BR*CR*CI - AI *BI - 2*AI *BI *CI 2 2 2 2 2 2 2 2 - 2*AI *BI *CR + 6*AI *BI *CI - 2*AI *BI*CR *CI 2 3 2 4 2 2 2 2 4 - 2*AI *BI*CI - AI *CR - 2*AI *CR *CI - AI *CI 4 3 2 2 + 2*AI*BR *CI + 4*AI*BR *CR*CI + 4*AI*BR *BI *CI 2 2 2 2 2 2 + 2*AI*BR *BI*CR - 2*AI*BR *BI*CI + 2*AI*BR *CR *CI 2 3 2 3 + 2*AI*BR *CI + 4*AI*BR*BI *CR*CI + 4*AI*BR*BI*CR 2 4 3 2 + 4*AI*BR*BI*CR*CI + 2*AI*BI *CI + 2*AI*BI *CR 3 2 2 2 2 3 - 2*AI*BI *CI - 2*AI*BI *CR *CI - 2*AI*BI *CI 4 2 2 4 4 2 + 2*AI*BI*CR + 4*AI*BI*CR *CI + 2*AI*BI*CI - BR *CR 4 2 3 3 3 2 2 2 2 - BR *CI - 2*BR *CR - 2*BR *CR*CI - 2*BR *BI *CR 2 2 2 2 2 2 3 2 4 - 2*BR *BI *CI + 2*BR *BI*CR *CI + 2*BR *BI*CI - BR *CR 2 2 2 2 4 2 3 2 2 - 2*BR *CR *CI - BR *CI - 2*BR*BI *CR - 2*BR*BI *CR*CI 4 2 4 2 3 2 3 3 2 4 - BI *CR - BI *CI + 2*BI *CR *CI + 2*BI *CI - BI *CR 2 2 2 2 4 - 2*BI *CR *CI - BI *CI ) > 0 COND clear x,x0,x1,x2,x3; %*********************************************************************** %***** ***** %***** T e s t Examples --- Module L I N B A N D ***** %***** ***** %*********************************************************************** % Example L.1 operator v; off echo; DIMENSION U(200),V(200),ACOF(200,3),ARHS(200),XL(200,3) DX=5.0E-2 X=1.0E-1 DO 25001 I=1,101 V(I)=X**2/2.0 X=X+DX 25001 CONTINUE IACOF=200 IARHS=200 N=1 AD(N)=1.0 AUCD(N)=0.0 ARHS(N)=V(1) N=N+1 ALCD(N)=1.0 AD(N)=-2.0 AUCD(N)=1.0 ARHS(N)=V(3)-(2.0*V(2))+V(1) DO 25002 K=3,99,1 N=N+1 ALCD(N)=1.0 AD(N)=-2.0 AUCD(N)=1.0 ARHS(N)=V(K-1)+V(K+1)-(2.0*V(K)) 25002 CONTINUE N=N+1 ALCD(N)=1.0 AD(N)=-2.0 AUCD(N)=1.0 ARHS(N)=V(101)-(2.0*V(100))+V(99) N=N+1 ALCD(N)=0.0 AD(N)=1.0 ARHS(N)=V(101) CALL DGTSL(N,ALCD,AD,AUCD,ARHS,IER) C N IS NUMBER OF EQUATIONS C ALCD,AD,AUCD,ARHS ARE ARRAYS OF DIMENSION AT LEAST (N) C IF (IER.NE.0) MATRIX ACOF IS ALGORITHMICALLY SINGULAR IF(IER.NE.0) CALL ERROUT N=1 U(1)=ARHS(N) N=N+1 U(2)=ARHS(N) DO 25003 K=3,99,1 N=N+1 U(K)=ARHS(N) 25003 CONTINUE N=N+1 U(100)=ARHS(N) N=N+1 U(101)=ARHS(N) AMER=0.0E-1 ARER=0.0E-1 DO 25004 I=1,101 AM=ABS(REAL(U(I)-V(I))) AR=AM/V(I) IF(AM.GT.AMER) AMER=AM IF(AR.GT.ARER) ARER=AR 25004 CONTINUE WRITE(*,100)AMER,ARER STOP 100 FORMAT(' MAX. ABS. ERROR = ',E12.2,' MAX. REL. ERROR = ',E12.2) END %*********************************************************************** % Example L.2 on nag; off echo; DIMENSION U(200),V(200),ACOF(200,3),ARHS(200),XL(200,3) DX=5.0E-2 X=1.0E-1 DO 25005 I=1,101 V(I)=X**2/2.0 X=X+DX 25005 CONTINUE IACOF=200 IARHS=200 N=1 AD(N)=1.0 AUCD(N+1)=0.0 ARHS(N)=V(1) N=N+1 ALCD(N)=1.0 AD(N)=-2.0 AUCD(N+1)=1.0 ARHS(N)=V(3)-(2.0*V(2))+V(1) DO 25006 K=3,99,1 N=N+1 ALCD(N)=1.0 AD(N)=-2.0 AUCD(N+1)=1.0 ARHS(N)=V(K-1)+V(K+1)-(2.0*V(K)) 25006 CONTINUE N=N+1 ALCD(N)=1.0 AD(N)=-2.0 AUCD(N+1)=1.0 ARHS(N)=V(101)-(2.0*V(100))+V(99) N=N+1 ALCD(N)=0.0 AD(N)=1.0 ARHS(N)=V(101) IER=0 CALL F01LEF(N,AD,0.,AUCD,ALCD,1.E-10,AU2CD,IN,IER) C N IS NUMBER OF EQUATIONS C ALCD,AD,AUCD,AU2CD,ARHS ARE ARRAYS OF DIMENSION AT LEAST (N) C IN IS INTEGER ARRAY OF DIMENSION AT LEAST (N) IF(IER.NE.0 .OR. IN(N).NE.0) CALL ERROUT CALL F04LEF(1,N,AD,AUCD,ALCD,AU2CD,IN,ARHS,0.,IER) IF(IER.NE.0) CALL ERROUT N=1 U(1)=ARHS(N) N=N+1 U(2)=ARHS(N) DO 25007 K=3,99,1 N=N+1 U(K)=ARHS(N) 25007 CONTINUE N=N+1 U(100)=ARHS(N) N=N+1 U(101)=ARHS(N) AMER=0.0E-1 ARER=0.0E-1 DO 25008 I=1,101 AM=ABS(REAL(U(I)-V(I))) AR=AM/V(I) IF(AM.GT.AMER) AMER=AM IF(AR.GT.ARER) ARER=AR 25008 CONTINUE WRITE(*,100)AMER,ARER STOP 100 FORMAT(' MAX. ABS. ERROR = ',E12.2,' MAX. REL. ERROR = ',E12.2) END %*********************************************************************** % Example L.3 on imsl; off echo,nag; DIMENSION U(200),V(200),ACOF(200,3),ARHS(200),XL(200,3) DX=5.0E-2 X=1.0E-1 DO 25009 I=1,101 V(I)=X**2/2.0 X=X+DX 25009 CONTINUE IACOF=200 IARHS=200 N=1 ACOF(N,1)=0.0 ACOF(N,2)=1.0 ACOF(N,3)=0.0 ARHS(N)=V(1) N=N+1 ACOF(N,1)=1.0 ACOF(N,2)=-2.0 ACOF(N,3)=1.0 ARHS(N)=V(3)-(2.0*V(2))+V(1) DO 25010 K=3,99,1 N=N+1 ACOF(N,1)=1.0 ACOF(N,2)=-2.0 ACOF(N,3)=1.0 ARHS(N)=V(K-1)+V(K+1)-(2.0*V(K)) 25010 CONTINUE N=N+1 ACOF(N,1)=1.0 ACOF(N,2)=-2.0 ACOF(N,3)=1.0 ARHS(N)=V(101)-(2.0*V(100))+V(99) N=N+1 ACOF(N,1)=0.0 ACOF(N,2)=1.0 ACOF(N,3)=0.0 ARHS(N)=V(101) CALL LEQT1B(ACOF,N,1,1,IACOF,ARHS,1,IARHS,0,XL,IER) C IACOF IS ACTUAL 1-ST DIMENSION OF THE ACOF ARRAY C IARHS IS ACTUAL 1-ST DIMENSION OF THE ARHS ARRAY C XL IS WORKING ARRAY WITH SIZE N*(NLC+1) C WHERE N IS NUMBER OF EQUATIONS NLC NUMBER OF LOWER C CODIAGONALS C IF IER=129( .NE.0) MATRIX ACOF IS ALGORITHMICALLY SINGULAR IF(IER.NE.0) CALL ERROUT N=1 U(1)=ARHS(N) N=N+1 U(2)=ARHS(N) DO 25011 K=3,99,1 N=N+1 U(K)=ARHS(N) 25011 CONTINUE N=N+1 U(100)=ARHS(N) N=N+1 U(101)=ARHS(N) AMER=0.0E-1 ARER=0.0E-1 DO 25012 I=1,101 AM=ABS(REAL(U(I)-V(I))) AR=AM/V(I) IF(AM.GT.AMER) AMER=AM IF(AR.GT.ARER) ARER=AR 25012 CONTINUE WRITE(*,100)AMER,ARER STOP 100 FORMAT(' MAX. ABS. ERROR = ',E12.2,' MAX. REL. ERROR = ',E12.2) END %*********************************************************************** % Example L.4 on essl; off echo,imsl; DIMENSION U(200),V(200),ACOF(200,3),ARHS(200),XL(200,3) DX=5.0E-2 X=1.0E-1 DO 25013 I=1,101 V(I)=X**2/2.0 X=X+DX 25013 CONTINUE IACOF=200 IARHS=200 N=1 AD(N)=1.0 AUCD(N)=0.0 ARHS(N)=V(1) N=N+1 ALCD(N)=1.0 AD(N)=-2.0 AUCD(N)=1.0 ARHS(N)=V(3)-(2.0*V(2))+V(1) DO 25014 K=3,99,1 N=N+1 ALCD(N)=1.0 AD(N)=-2.0 AUCD(N)=1.0 ARHS(N)=V(K-1)+V(K+1)-(2.0*V(K)) 25014 CONTINUE N=N+1 ALCD(N)=1.0 AD(N)=-2.0 AUCD(N)=1.0 ARHS(N)=V(101)-(2.0*V(100))+V(99) N=N+1 ALCD(N)=0.0 AD(N)=1.0 ARHS(N)=V(101) CALL DGTF(N,ALCD,AD,AUCD,AF,IPVT) C N IS NUMBER OF EQUATIONS C ALCD,AD,AUCD,AF,ARHS ARE ARRAYS OF DIMENSION AT LEAST (N) C THESE ARRAYS ARE DOUBLE PRECISION TYPE C FOR SINGLE PRECISION CHANGE DGTF TO SGTF AND DGTS TO SGTS C IPVT IS INTEGER ARRAY OF DIMENSION AT LEAST (N+3)/8 CALL DGTS(N,ALCD,AD,AUCD,AF,IPVT,ARHS) N=1 U(1)=ARHS(N) N=N+1 U(2)=ARHS(N) DO 25015 K=3,99,1 N=N+1 U(K)=ARHS(N) 25015 CONTINUE N=N+1 U(100)=ARHS(N) N=N+1 U(101)=ARHS(N) AMER=0.0E-1 ARER=0.0E-1 DO 25016 I=1,101 AM=ABS(REAL(U(I)-V(I))) AR=AM/V(I) IF(AM.GT.AMER) AMER=AM IF(AR.GT.ARER) ARER=AR 25016 CONTINUE WRITE(*,100)AMER,ARER STOP 100 FORMAT(' MAX. ABS. ERROR = ',E12.2,' MAX. REL. ERROR = ',E12.2) END off essl; %*********************************************************************** %***** ***** %***** T e s t Complex Examples --- More Modules ***** %***** ***** %*********************************************************************** % Example M.1 off exp; coordinates t,x into n,j; grid uniform,x,t; dependence v(t,x),w(t,x); %isgrid u(x..half),v(x..half); %grideq u(x..one),v(x..one); iim aa, v, diff(v,t)=c*diff(w,x), w, diff(w,t)=c*diff(v,x); ***************************** ***** Program ***** IIMET Ver 1.1 17/05/91 ***************************** Partial Differential Equations ============================== DIFF(V,T) - DIFF(W,X)*C = 0 - DIFF(V,X)*C + DIFF(W,T) = 0 0 interpolations are needed in T coordinate Equation for V variable is integrated in HALF grid point Equation for W variable is integrated in HALF grid point 0 interpolations are needed in X coordinate Equation for V variable is integrated in HALF grid point Equation for W variable is integrated in ONE grid point Equations after Discretization Using IIM : ========================================== ( - (W(N + 1,J + 1) - W(N + 1,J) + W(N,J + 1) - W(N,J))*C*HT 2*J + 1 2*J + 1 + 2*V(N + 1,---------)*HX - 2*V(N,---------)*HX)/(2*HX*HT) = 0 2 2 2*J - 1 2*J + 1 2*J - 1 ((V(N + 1,---------) - V(N + 1,---------) + V(N,---------) 2 2 2 2*J + 1 - V(N,---------))*C*HT + 2*W(N + 1,J)*HX - 2*W(N,J)*HX)/(2*HX*HT) 2 = 0 on exp; center t=1/2; functions v,w; approx( aa(0,0)=aa(0,1)); Difference scheme approximates differential equation DF(V,T) - DF(W,X)*C=0 with orders of approximation: 2 HT HX center x=1/2; approx( aa(1,0)=aa(1,1)); Difference scheme approximates differential equation - DF(V,X)*C + DF(W,T)=0 with orders of approximation: 2 HT HX let cos ax**2=1-sin ax**2; unfunc v,w; matrix a(1,2),b(2,2),bt(2,2); a(1,1):=aa(0,0); 2*J + 1 2*J + 1 A(1,1) := (2*V(N + 1,---------)*HX - 2*V(N,---------)*HX 2 2 - W(N + 1,J + 1)*C*HT + W(N + 1,J)*C*HT - W(N,J + 1)*C*HT + W(N,J)*C*HT)/(2*HX*HT) a(1,2):=aa(1,0); 2*J - 1 2*J + 1 A(1,2) := (V(N + 1,---------)*C*HT - V(N + 1,---------)*C*HT 2 2 2*J - 1 2*J + 1 + V(N,---------)*C*HT - V(N,---------)*C*HT 2 2 + 2*W(N + 1,J)*HX - 2*W(N,J)*HX)/(2*HX*HT) off prfourmat; b:=ampmat a; KX*HX AX := ------- 2 [ 2 2 2 2 ] [ - SIN(AX) *C *HT + HX 2*SIN(AX)*C*I*HX*HT ] [-------------------------- ----------------------- ] [ 2 2 2 2 2 2 2 2 ] [ SIN(AX) *C *HT + HX SIN(AX) *C *HT + HX ] B := [ ] [ 2 2 2 2 ] [ 2*SIN(AX)*C*I*HX*HT - SIN(AX) *C *HT + HX ] [ ----------------------- --------------------------] [ 2 2 2 2 2 2 2 2 ] [ SIN(AX) *C *HT + HX SIN(AX) *C *HT + HX ] clear a,aa; factor lam; pol:=charpol b; 2 4 4 4 2 2 2 2 4 POL := (LAM *(SIN(AX) *C *HT + 2*SIN(AX) *C *HX *HT + HX ) 4 4 4 4 4 4 4 + 2*LAM*(SIN(AX) *C *HT - HX ) + SIN(AX) *C *HT 2 2 2 2 4 4 4 4 + 2*SIN(AX) *C *HX *HT + HX )/(SIN(AX) *C *HT 2 2 2 2 4 + 2*SIN(AX) *C *HX *HT + HX ) pol:=troot1 pol; 2 2 2 4*SIN(AX) *C *HT If ----------------------- = 0 and 0 2 2 2 2 SIN(AX) *C *HT + HX = 0 , a root of the polynomial is equal to 1. 2 4 4 4 2 2 2 2 4 POL := (LAM *(SIN(AX) *C *HT + 2*SIN(AX) *C *HX *HT + HX ) 4 4 4 4 4 4 4 + 2*LAM*(SIN(AX) *C *HT - HX ) + SIN(AX) *C *HT 2 2 2 2 4 4 4 4 + 2*SIN(AX) *C *HX *HT + HX )/(SIN(AX) *C *HT 2 2 2 2 4 + 2*SIN(AX) *C *HX *HT + HX ) pol:=hurw num pol; 2 2 2 2 2 2 2 2 POL := 4*LAM *SIN(AX) *C *HT *(SIN(AX) *C *HT + HX ) 2 2 2 2 2 + 4*HX *(SIN(AX) *C *HT + HX ) hurwitzp pol; Necessary and sufficient conditions are: 2 HX ----------------- > 0 2 2 2 SIN(AX) *C *HT NO bt:=tcon b; [ 2 2 2 2 ] [ - SIN(AX) *C *HT + HX - 2*SIN(AX)*C*I*HX*HT ] [-------------------------- ------------------------ ] [ 2 2 2 2 2 2 2 2 ] [ SIN(AX) *C *HT + HX SIN(AX) *C *HT + HX ] BT := [ ] [ 2 2 2 2 ] [ - 2*SIN(AX)*C*I*HX*HT - SIN(AX) *C *HT + HX ] [ ------------------------ --------------------------] [ 2 2 2 2 2 2 2 2 ] [ SIN(AX) *C *HT + HX SIN(AX) *C *HT + HX ] bt*b; [1 0] [ ] [0 1] bt*b-b*bt; [0 0] [ ] [0 0] clear aa,a,b,bt; %*********************************************************************** % Example M.2 : Richtmyer, Morton: Difference methods for initial value % problems, &10.2. p.261 coordinates t,x into n,j; grid uniform,t,x; let cos ax**2=1-sin ax**2; unfunc v,w; matrix a(1,2),b(2,2),bt(2,2); a(1,1):=(v(n+1)-v(n))/ht-c*(w(j+1/2)-w(j-1/2))/hx$ a(1,2):=(w(n+1,j-1/2)-w(n,j-1/2))/ht-c*(v(n+1,j)-v(n+1,j-1))/hx$ off prfourmat; b:=ampmat a; KX*HX AX := ------- 2 [ 2*SIN(AX)*C*I*HT ] [ 1 ------------------ ] [ HX ] [ ] B := [ 2 2 2 2 ] [ 2*SIN(AX)*C*I*HT - 4*SIN(AX) *C *HT + HX ] [------------------ ----------------------------] [ HX 2 ] [ HX ] clear a; factor lam; pol:=charpol b; 2 2 2 2 2 2 2 LAM *HX + 2*LAM*(2*SIN(AX) *C *HT - HX ) + HX POL := -------------------------------------------------- 2 HX pol:=hurw num pol; 2 2 2 2 2 2 2 2 POL := 4*LAM *SIN(AX) *C *HT + 4*( - SIN(AX) *C *HT + HX ) hurwitzp pol; Necessary and sufficient conditions are: 2 2 2 2 - SIN(AX) *C *HT + HX -------------------------- > 0 2 2 2 SIN(AX) *C *HT NO bt:=tcon b; [ - 2*SIN(AX)*C*I*HT ] [ 1 --------------------- ] [ HX ] [ ] BT := [ 2 2 2 2 ] [ - 2*SIN(AX)*C*I*HT - 4*SIN(AX) *C *HT + HX ] [--------------------- ----------------------------] [ HX 2 ] [ HX ] bt*b; 2 2 2 2 3 3 3 4*SIN(AX) *C *HT + HX 8*SIN(AX) *C *I*HT MAT((-------------------------,---------------------), 2 3 HX HX 3 3 3 - 8*SIN(AX) *C *I*HT (------------------------, 3 HX 4 4 4 2 2 2 2 4 16*SIN(AX) *C *HT - 4*SIN(AX) *C *HX *HT + HX --------------------------------------------------)) 4 HX bt*b-b*bt; [ 3 3 3 ] [ 16*SIN(AX) *C *I*HT ] [ 0 ----------------------] [ 3 ] [ HX ] [ ] [ 3 3 3 ] [ - 16*SIN(AX) *C *I*HT ] [------------------------- 0 ] [ 3 ] [ HX ] clear a,b,bt; %*********************************************************************** % Example M.3: Mazurik: Algoritmy resenia zadaci..., preprint no.24-85, % AN USSR SO, Inst. teor. i prikl. mechaniky, p.34 operator v1,v2; matrix a(1,3),b(3,3),bt(3,3); a(1,1):=(p(n+1)-p(n))/ht+c/2*((-p(m-1)+2*p(m)-p(m+1))/hx + (v1(m+1)-v1(m-1))/hx - (p(k-1)-2*p(k)+p(k+1))/hy + (v2(k+1)-v2(k-1))/hy)$ a(1,2):=(v1(n+1)-v1(n))/ht+c/2*((p(m+1)-p(m-1))/hx - (v1(m-1)-2*v1(m)+v1(m+1))/hx)$ a(1,3):=(v2(n+1)-v2(n))/ht + c/2*((p(k+1)-p(k-1))/hy - (v2(k-1)-2*v2(k)+v2(k+1))/hy)$ coordinates t,x,y into n,m,k; functions p,v1,v2; for k:=1:3 do approx(a(1,k)=0); Difference scheme approximates differential equation DF(P,T) + DF(V1,X)*C + DF(V2,Y)*C=0 with orders of approximation: 2 HT HX HY Difference scheme approximates differential equation DF(P,X)*C + DF(V1,T)=0 with orders of approximation: 2 HT HX 1 Difference scheme approximates differential equation DF(P,Y)*C + DF(V2,T)=0 with orders of approximation: 2 HT 1 HY grid uniform,t,x,y; unfunc p,v1,v2; hy:=hx; HY := HX off prfourmat; b:=ampmat a; AX := KX*HX AY := KY*HY COS(AX)*C*HT + COS(AY)*C*HT - 2*C*HT + HX B := MAT((-------------------------------------------, HX - SIN(AX)*C*I*HT - SIN(AY)*C*I*HT -------------------,-------------------), HX HX - SIN(AX)*C*I*HT COS(AX)*C*HT - C*HT + HX (-------------------,--------------------------,0), HX HX - SIN(AY)*C*I*HT COS(AY)*C*HT - C*HT + HX (-------------------,0,--------------------------)) HX HX pol:=charpol b; 3 3 2 2 POL := (LAM *HX + LAM *HX *( - 2*COS(AX)*C*HT - 2*COS(AY)*C*HT + 4*C*HT - 3*HX) + LAM 2 2 2 2 *HX*(3*COS(AX)*COS(AY)*C *HT - 5*COS(AX)*C *HT 2 2 + 4*COS(AX)*C*HX*HT - 5*COS(AY)*C *HT 2 2 2 + 4*COS(AY)*C*HX*HT + 7*C *HT - 8*C*HX*HT + 3*HX ) 3 3 2 2 + 4*COS(AX)*COS(AY)*C *HT - 3*COS(AX)*COS(AY)*C *HX*HT 3 3 2 2 - 4*COS(AX)*C *HT + 5*COS(AX)*C *HX*HT 2 3 3 - 2*COS(AX)*C*HX *HT - 4*COS(AY)*C *HT 2 2 2 3 3 + 5*COS(AY)*C *HX*HT - 2*COS(AY)*C*HX *HT + 4*C *HT 2 2 2 3 3 - 7*C *HX*HT + 4*C*HX *HT - HX )/HX let cos ax=cos ax2**2-sin ax2**2, cos ay=cos ay2**2-sin ay2**2, sin ax=2*sin ax2*cos ax2, sin ay=2*sin ay2*cos ay2, cos ax2**2=1-sin ax2**2, cos ay2**2=1-sin ay2**2, sin ax2=s1, sin ay2=s2, hx=c*ht/cap; factor lam; order s1,s2; pol:=troot1 pol; 2 2 3 If 16*S1 *S2 *CAP = 0 and 0 = 0 , a root of the polynomial is equal to 1. 3 2 2 2 2 2 2 POL := LAM + LAM *(4*S1 *CAP + 4*S2 *CAP - 3) + LAM*(12*S1 *S2 *CAP 2 2 2 2 2 2 + 4*S1 *CAP - 8*S1 *CAP + 4*S2 *CAP - 8*S2 *CAP + 3) 2 2 3 2 2 2 2 2 2 + 16*S1 *S2 *CAP - 12*S1 *S2 *CAP - 4*S1 *CAP + 4*S1 *CAP 2 2 2 - 4*S2 *CAP + 4*S2 *CAP - 1 clear cos ax,cos ay,sin ax,sin ay,cos ax2**2,cos ay2**2,sin ax2,sin ay2, hx,hy; pol:=hurw num pol; 3 2 2 3 POL := 16*LAM *S1 *S2 *CAP 2 2 2 2 2 2 2 2 + 8*LAM *CAP *( - 6*S1 *S2 *CAP + 3*S1 *S2 + S1 + S2 ) + 16 2 2 2 2 2 2 2 *LAM*CAP*(3*S1 *S2 *CAP - 3*S1 *S2 *CAP - S1 *CAP + S1 2 2 2 2 3 - S2 *CAP + S2 ) + 8*( - 2*S1 *S2 *CAP 2 2 2 2 2 2 2 2 + 3*S1 *S2 *CAP + S1 *CAP - 2*S1 *CAP + S2 *CAP 2 - 2*S2 *CAP + 1) hurwitzp pol; If we denote: 2 2 3 (1) 16*S1 *S2 *CAP > 0 2 2 2 2 2 2 2 (2) 8*CAP *( - 6*S1 *S2 *CAP + 3*S1 *S2 + S1 + S2 ) > 0 (3) 16*CAP 2 2 2 2 2 2 2 2 2 *(3*S1 *S2 *CAP - 3*S1 *S2 *CAP - S1 *CAP + S1 - S2 *CAP + S2 ) > 0 2 2 3 2 2 2 2 2 2 (4) 8*( - 2*S1 *S2 *CAP + 3*S1 *S2 *CAP + S1 *CAP - 2*S1 *CAP 2 2 2 + S2 *CAP - 2*S2 *CAP + 1) > 0 2 2 2 2 2 2 2 (5) 8*CAP *( - 6*S1 *S2 *CAP + 3*S1 *S2 + S1 + S2 ) > 0 3 4 4 3 4 4 2 4 4 (6) 128*CAP *( - 16*S1 *S2 *CAP + 24*S1 *S2 *CAP - 9*S1 *S2 *CAP 4 2 2 4 2 4 2 + 8*S1 *S2 *CAP - 10*S1 *S2 *CAP + 3*S1 *S2 4 4 2 4 2 2 4 - S1 *CAP + S1 + 8*S1 *S2 *CAP - 10*S1 *S2 *CAP 2 4 2 2 2 2 4 + 3*S1 *S2 - 2*S1 *S2 *CAP + S1 *S2 - S2 *CAP 4 + S2 ) > 0 3 6 6 6 6 6 5 6 6 4 (7) 1024*CAP *(32*S1 *S2 *CAP - 96*S1 *S2 *CAP + 90*S1 *S2 *CAP 6 6 3 6 4 5 - 27*S1 *S2 *CAP - 32*S1 *S2 *CAP 6 4 4 6 4 3 + 100*S1 *S2 *CAP - 93*S1 *S2 *CAP 6 4 2 6 2 4 + 27*S1 *S2 *CAP + 10*S1 *S2 *CAP 6 2 3 6 2 2 6 2 - 31*S1 *S2 *CAP + 26*S1 *S2 *CAP - 6*S1 *S2 *CAP 6 3 6 2 6 - S1 *CAP + 3*S1 *CAP - 2*S1 *CAP 4 6 5 4 6 4 - 32*S1 *S2 *CAP + 100*S1 *S2 *CAP 4 6 3 4 6 2 - 93*S1 *S2 *CAP + 27*S1 *S2 *CAP 4 4 4 4 4 3 + 20*S1 *S2 *CAP - 76*S1 *S2 *CAP 4 4 2 4 4 4 2 3 + 73*S1 *S2 *CAP - 21*S1 *S2 *CAP - 3*S1 *S2 *CAP 4 2 2 4 2 4 2 + 16*S1 *S2 *CAP - 14*S1 *S2 *CAP + 3*S1 *S2 4 4 2 6 4 2 6 3 - S1 *CAP + S1 + 10*S1 *S2 *CAP - 31*S1 *S2 *CAP 2 6 2 2 6 2 4 3 + 26*S1 *S2 *CAP - 6*S1 *S2 *CAP - 3*S1 *S2 *CAP 2 4 2 2 4 2 4 + 16*S1 *S2 *CAP - 14*S1 *S2 *CAP + 3*S1 *S2 2 2 2 2 6 3 6 2 - 2*S1 *S2 *CAP + S1 *S2 - S2 *CAP + 3*S2 *CAP 6 4 4 - 2*S2 *CAP - S2 *CAP + S2 ) > 0 (C1) (1) AND (2) AND (4) (C2) (1) AND (3) AND (4) (D1) (5) AND (7) (D2) (6) Necessary and sufficient conditions are: ( (C1) OR (C2) ) AND ( (D1) OR (D2) ) COND bt:=tcon b; COS(AX)*C*HT + COS(AY)*C*HT - 2*C*HT + HX BT := MAT((-------------------------------------------, HX SIN(AX)*C*I*HT SIN(AY)*C*I*HT ----------------,----------------), HX HX SIN(AX)*C*I*HT COS(AX)*C*HT - C*HT + HX (----------------,--------------------------,0), HX HX SIN(AY)*C*I*HT COS(AY)*C*HT - C*HT + HX (----------------,0,--------------------------)) HX HX bt*b; 2 2 2 2 MAT(((2*COS(AX)*COS(AY)*C *HT - 4*COS(AX)*C *HT + 2*COS(AX)*C*HX*HT 2 2 2 2 - 4*COS(AY)*C *HT + 2*COS(AY)*C*HX*HT + 6*C *HT - 4*C*HX*HT 2 2 2 2 SIN(AX)*C *I*HT *( - COS(AY) + 1) + HX )/HX ,-----------------------------------, 2 HX 2 2 SIN(AY)*C *I*HT *( - COS(AX) + 1) -----------------------------------), 2 HX 2 2 SIN(AX)*C *I*HT *(COS(AY) - 1) 2 2 (--------------------------------,( - 2*COS(AX)*C *HT 2 HX 2 2 2 2 + 2*COS(AX)*C*HX*HT + 2*C *HT - 2*C*HX*HT + HX )/HX , 2 2 SIN(AX)*SIN(AY)*C *HT ------------------------), 2 HX 2 2 2 2 SIN(AY)*C *I*HT *(COS(AX) - 1) SIN(AX)*SIN(AY)*C *HT (--------------------------------,------------------------,( 2 2 HX HX 2 2 2 2 - 2*COS(AY)*C *HT + 2*COS(AY)*C*HX*HT + 2*C *HT 2 2 - 2*C*HX*HT + HX )/HX )) bt*b-b*bt; 2 2 2*SIN(AX)*C *I*HT *( - COS(AY) + 1) MAT((0,-------------------------------------, 2 HX 2 2 2*SIN(AY)*C *I*HT *( - COS(AX) + 1) -------------------------------------), 2 HX 2 2 2*SIN(AX)*C *I*HT *(COS(AY) - 1) (----------------------------------,0,0), 2 HX 2 2 2*SIN(AY)*C *I*HT *(COS(AX) - 1) (----------------------------------,0,0)) 2 HX clear a,b,bt,pol; %*********************************************************************** end; Quitting