@@ -1,3803 +1,3803 @@ -REDUCE 3.6, 15-Jul-95, patched to 6 Mar 96 ... - - -%*********************************************************************** -%***** ***** -%***** Package F I D E - Test Examples Ver. 1.1.2 May 29,1995 ***** -%***** ***** -%*********************************************************************** -%*********************************************************************** -%***** ***** -%***** 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; - - -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) cos(th)*u(2) + 2*sin(th)*u(1) -a2:=------------- + ------------- + df(u(1),r) + ------------------------------- - sin(th)*r r sin(th)*r - - -a3:=curl v; - - df(v(3),th) - df(v(2),fi) cos(th)*v(3) -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) df(v(1),fi,2) -a4:= ( ------------------ + ------------------ + --------------- + df(v(1),r,2) - 2 2 2 2 - sin(th)*r r sin(th) *r - - 2*df(v(1),r) df(v(1),th,2) df(v(1),th)*cos(th) - + -------------- + --------------- + --------------------- - r 2 2 - r sin(th)*r - - 2*(cos(th)*v(2) + sin(th)*v(1)) - - --------------------------------- , - 2 - sin(th)*r - - - 2*df(v(3),fi)*cos(th) df(v(2),fi,2) - -------------------------- + --------------- + df(v(2),r,2) - 2 2 2 2 - sin(th) *r 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 - - 2*df(v(1),th) - v(2) - + --------------- + ------------- , - 2 2 2 - r sin(th) *r - - df(v(3),fi,2) 2*df(v(3),r) df(v(3),th,2) - --------------- + df(v(3),r,2) + -------------- + --------------- - 2 2 r 2 - sin(th) *r r - - df(v(3),th)*cos(th) 2*df(v(2),fi)*cos(th) 2*df(v(1),fi) - + --------------------- + ----------------------- + --------------- - 2 2 2 2 - sin(th)*r sin(th) *r sin(th)*r - - - v(3) - + ------------- ) - 2 2 - sin(th) *r - - -a3:=2*a3+a4; - - - 2*df(v(3),fi) 2*df(v(3),th) - 2*df(v(2),fi) -a3:= ( ------------------ + --------------- + ------------------ - 2 r sin(th)*r - sin(th)*r - - - 2*df(v(2),th) df(v(1),fi,2) 2*df(v(1),r) - + ------------------ + --------------- + df(v(1),r,2) + -------------- - 2 2 2 r - r sin(th) *r - - df(v(1),th,2) df(v(1),th)*cos(th) - + --------------- + --------------------- - 2 2 - r sin(th)*r - - 2*(cos(th)*v(3)*r - cos(th)*v(2) - sin(th)*v(1)) - + -------------------------------------------------- , - 2 - sin(th)*r - - - 2*df(v(3),fi)*cos(th) df(v(2),fi,2) - -------------------------- - 2*df(v(3),r) + --------------- - 2 2 2 2 - sin(th) *r sin(th) *r - - 2*df(v(2),r) df(v(2),th,2) - + df(v(2),r,2) + -------------- + --------------- - r 2 - r - - df(v(2),th)*cos(th) 2*df(v(1),fi) 2*df(v(1),th) - + --------------------- + --------------- + --------------- - 2 sin(th)*r 2 - sin(th)*r r - - 2 - - 2*sin(th) *v(3)*r - v(2) - + ----------------------------- , - 2 2 - sin(th) *r - - df(v(3),fi,2) 2*df(v(3),r) df(v(3),th,2) - --------------- + df(v(3),r,2) + -------------- + --------------- - 2 2 r 2 - sin(th) *r r - - df(v(3),th)*cos(th) 2*df(v(2),fi)*cos(th) - + --------------------- + ----------------------- + 2*df(v(2),r) - 2 2 2 - sin(th)*r sin(th) *r - - 2 - 2*df(v(1),fi) - 2*df(v(1),th) 2*sin(th) *v(2)*r - v(3) - + --------------- + ------------------ + -------------------------- ) - 2 r 2 2 - sin(th)*r sin(th) *r - - -a5:=lapl f; - - df(f,fi,2) 2*df(f,r) df(f,th,2) df(f,th)*cos(th) -a5:=------------- + df(f,r,2) + ----------- + ------------ + ------------------ - 2 2 r 2 2 - sin(th) *r r 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 - - cos(th)*w(2,1) - sin(th)*w(3,3) - sin(th)*w(2,2) + 2*sin(th)*w(1,1) - + --------------------------------------------------------------------- - 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 - - - cos(th)*w(3,3) + cos(th)*w(2,2) + sin(th)*w(2,1) + 2*sin(th)*w(1,2) - ------------------------------------------------------------------------ - 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 - - cos(th)*w(3,2) + cos(th)*w(2,3) + sin(th)*w(3,1) + 2*sin(th)*w(1,3) - + --------------------------------------------------------------------- - 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),fi)*u(3) df(v(1),th)*u(2) - ( ------------------ + df(v(1),r)*u(1) + ------------------ - sin(th)*r r - - u(3)*v(3) + u(2)*v(2) - - ----------------------- , - r - - df(v(2),fi)*u(3) df(v(2),th)*u(2) - ------------------ + df(v(2),r)*u(1) + ------------------ - sin(th)*r r - - - cos(th)*u(3)*v(3) + sin(th)*u(2)*v(1) - + ------------------------------------------ , - sin(th)*r - - df(v(3),fi)*u(3) df(v(3),th)*u(2) - ------------------ + df(v(3),r)*u(1) + ------------------ - sin(th)*r r - - u(3)*(cos(th)*v(2) + sin(th)*v(1)) - + ------------------------------------ ) - 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; - - - -%*********************************************************************** -%***** ***** -%***** 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.2 -***************************** - - Partial Differential Equations - ============================== - -diff(eta,t) - diff(v,x)*eta = 0 - - diff(p,x)*eta ---------------- + diff(v,t) = 0 - ro - - diff(v,x)*eta*p -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 : - ========================================== - -(4*(eta(j,m + 1) - eta(j,m) - eta(j + 1,m) + eta(j + 1,m + 1))*hx - ( - - (eta(j + 1,m + 1) + eta(j,m + 1))*(v(j + 1,m + 1) - v(j,m + 1)) - - + (eta(j + 1,m) + eta(j,m))*(v(j + 1,m) - v(j,m)))*(ht(m + 1) + ht(m)))/(4 - - *(ht(m + 1) + ht(m))*hx) = 0 - - -(4*(v(j,m + 1) - v(j,m) - v(j + 1,m) + v(j + 1,m + 1))*hx*ro + ( - - (eta(j + 1,m + 1) + eta(j,m + 1))*(p(j + 1,m + 1) - p(j,m + 1)) - - + (eta(j + 1,m) + eta(j,m))*(p(j + 1,m) - p(j,m)))*(ht(m + 1) + ht(m)))/(4 - - *(ht(m + 1) + ht(m))*hx*ro) = 0 - - -(4*(eps(j,m + 1) - eps(j,m) - eps(j + 1,m) + eps(j + 1,m + 1))*hx*ro + ( - - (eta(j + 1,m + 1)*p(j + 1,m + 1) + eta(j,m + 1)*p(j,m + 1)) - - *(v(j + 1,m + 1) - v(j,m + 1)) - - + (eta(j + 1,m)*p(j + 1,m) + eta(j,m)*p(j,m))*(v(j + 1,m) - v(j,m))) - - *(ht(m + 1) + ht(m)))/(4*(ht(m + 1) + ht(m))*hx*ro) = 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.2 -***************************** - - 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 : - ========================================== - - 2*j + 1 2*j + 1 2*j + 3 - - ((2*(v(---------,m + 1) + v(---------,m)) - v(---------,m) - 2 2 2 - - 2*j + 3 2*j - 1 2*j - 1 - - v(---------,m + 1) - v(---------,m) - v(---------,m + 1)) - 2 2 2 - - 2*j + 1 - *cos(x(---------))*ht + ( - 2 - - (u(j,m + 1) + u(j,m) - u(j + 1,m) - u(j + 1,m + 1))*ht - - 2 - - (u(j,m + 1) - u(j,m) - u(j + 1,m) + u(j + 1,m + 1))*hx)*hx)/(2*ht*hx ) - - = 0 - - -( - ((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 2 - - + (u(j + 1,m + 1) + u(j + 1,m))*sin(x(j + 1))*ht)/(2*ht*hx) = 0 - - - -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.2 -***************************** - - Partial Differential Equations - ============================== - - diff(ur,x,2) 2 2 - - diff(ui,t) + -------------- = - ur*(ui + ur ) - 2 - - diff(ui,x,2) 2 2 --------------- + diff(ur,t) = - ui*(ui + ur ) - 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,j + 1) - 2*ur(m,j) + ur(m,j - 1) - 2*ur(m + 1,j) + ur(m + 1,j + 1) - - 2 2 - + ur(m + 1,j - 1))*ht - 4*(ui(m + 1,j) - ui(m,j))*hx )/(4*ht*hx ) = - - 3 3 2 2 - ur(m + 1,j) + ur(m,j) + ui(m,j) *ur(m,j) + ui(m + 1,j) *ur(m + 1,j) - - ----------------------------------------------------------------------- - 2 - - -((ui(m,j + 1) - 2*ui(m,j) + ui(m,j - 1) - 2*ui(m + 1,j) + ui(m + 1,j + 1) - - 2 2 - + ui(m + 1,j - 1))*ht + 4*(ur(m + 1,j) - ur(m,j))*hx )/(4*ht*hx ) = - - 2 2 3 2 - (ui(m,j) + ur(m,j) )*ui(m,j) + ui(m + 1,j) + ui(m + 1,j)*ur(m + 1,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; - - -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.2 -***************************** - - 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 5*diff(u1,x)*p -------------- + ---------------- + ---------------- + ---------------- - 2 2 2 2 - - 5*diff(u2,y)*p - + ---------------- = 0 - 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 *(n(j + 1,i + 1,m + 1)*u1(j + 1,i + 1,m + 1)*hy - - + n(j + 1,i + 1,m + 1)*u2(j + 1,i + 1,m + 1)*hx - - + n(j + 1,i + 1,m)*u1(j + 1,i + 1,m)*hy - - + n(j + 1,i + 1,m)*u2(j + 1,i + 1,m)*hx - - + n(j + 1,i,m + 1)*u1(j + 1,i,m + 1)*hy - - - n(j + 1,i,m + 1)*u2(j + 1,i,m + 1)*hx - - + n(j + 1,i,m)*u1(j + 1,i,m)*hy - n(j + 1,i,m)*u2(j + 1,i,m)*hx - - - n(j,i + 1,m + 1)*u1(j,i + 1,m + 1)*hy - - + n(j,i + 1,m + 1)*u2(j,i + 1,m + 1)*hx - - - n(j,i + 1,m)*u1(j,i + 1,m)*hy + n(j,i + 1,m)*u2(j,i + 1,m)*hx - - - n(j,i,m + 1)*u1(j,i,m + 1)*hy - n(j,i,m + 1)*u2(j,i,m + 1)*hx - - -1 - - n(j,i,m)*u1(j,i,m)*hy - n(j,i,m)*u2(j,i,m)*hx))/4 + (ht *( - - n(j,i,m + 1) - n(j,i,m) - n(j,i + 1,m) + n(j,i + 1,m + 1) - n(j + 1,i,m) - - + n(j + 1,i,m + 1) - n(j + 1,i + 1,m) + n(j + 1,i + 1,m + 1)))/4 = 0 - - - -1 -(hx *((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)) + - - (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 + 1,m + 1)*u1(j + 1,i + 1,m + 1) - - + n(j,i + 1,m + 1)*u1(j,i + 1,m + 1)) - - -1 -1 -1 - *(u1(j + 1,i + 1,m + 1) - u1(j,i + 1,m + 1)))*m)/8 + (hy *hx *ht *((( - - (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)) + - - (n(j + 1,i,m + 1) + n(j + 1,i,m)) - - *(u1(j + 1,i,m + 1) - u1(j + 1,i,m)) + - - (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)))*hx*m + 2*( - - p(j + 1,i,m + 1) + p(j + 1,i,m) + p(j + 1,i + 1,m) - - + p(j + 1,i + 1,m + 1) - - - (p(j,i,m + 1) + p(j,i,m) + p(j,i + 1,m) + p(j,i + 1,m + 1)))*ht) - - *hy + ((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)) - - *(u1(j,i + 1,m) - u1(j,i,m)) + - - (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 + 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)))*ht*hx*m))/8 = 0 - - - -1 -1 -(hy *hx *(((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)) - - *(u2(j + 1,i,m) - u2(j,i,m)) + - - (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 + 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)))*hy + ( - - (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)) - - *(u2(j,i + 1,m) - u2(j,i,m)) + - - (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 + 1,i + 1,m + 1)*u2(j + 1,i + 1,m + 1) - - + n(j + 1,i,m + 1)*u2(j + 1,i,m + 1)) - - -1 - *(u2(j + 1,i + 1,m + 1) - u2(j + 1,i,m + 1)))*hx)*m)/8 + ( - hy - - -1 - *ht *(2*(p(j,i,m + 1) + p(j,i,m) - p(j,i + 1,m) - p(j,i + 1,m + 1) - - + p(j + 1,i,m) + p(j + 1,i,m + 1) - p(j + 1,i + 1,m) - - - p(j + 1,i + 1,m + 1))*ht - ((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)) + - - (n(j + 1,i,m + 1) + n(j + 1,i,m)) - - *(u2(j + 1,i,m + 1) - u2(j + 1,i,m)) + - - (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)))*hy*m))/8 = 0 - - - -1 -1 -(hy *hx *(3*((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)) + - - (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 + 1,m + 1) - p(j,i + 1,m + 1)) - - *(u1(j + 1,i + 1,m + 1) + u1(j,i + 1,m + 1)))*hy + 2*( - - 4*p(j + 1,i + 1,m + 1)*u2(j + 1,i + 1,m + 1) - - - p(j + 1,i + 1,m + 1)*u2(j + 1,i,m + 1) - - + 4*p(j + 1,i + 1,m)*u2(j + 1,i + 1,m) - - - p(j + 1,i + 1,m)*u2(j + 1,i,m) - - + p(j + 1,i,m + 1)*u2(j + 1,i + 1,m + 1) - - - 4*p(j + 1,i,m + 1)*u2(j + 1,i,m + 1) - - + p(j + 1,i,m)*u2(j + 1,i + 1,m) - 4*p(j + 1,i,m)*u2(j + 1,i,m) - - + 4*p(j,i + 1,m + 1)*u2(j,i + 1,m + 1) - - - p(j,i + 1,m + 1)*u2(j,i,m + 1) + 4*p(j,i + 1,m)*u2(j,i + 1,m) - - - p(j,i + 1,m)*u2(j,i,m) + p(j,i,m + 1)*u2(j,i + 1,m + 1) - - - 4*p(j,i,m + 1)*u2(j,i,m + 1) + p(j,i,m)*u2(j,i + 1,m) - - - 4*p(j,i,m)*u2(j,i,m))*hx + 5*( - - (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)) + - - (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 + 1,m + 1) + p(j,i + 1,m + 1)) - - -1 - *(u1(j + 1,i + 1,m + 1) - u1(j,i + 1,m + 1)))*hy))/16 + (3*ht *( - - p(j,i,m + 1) - p(j,i,m) - p(j,i + 1,m) + p(j,i + 1,m + 1) - p(j + 1,i,m) - - + p(j + 1,i,m + 1) - p(j + 1,i + 1,m) + p(j + 1,i + 1,m + 1)))/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.2 -***************************** - - 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 - - 3*diff(tt,t)*k*n 3*diff(tt,x)*k*n*u -diff(q,x) + ------------------ + -------------------- - 2 2 - - diff(u,x)*(2*k*n*tt + p) - + -------------------------- = 0 - 2 - - 2*diff(q,x) -diff(p,t) + diff(p,x)*u + ------------- + diff(u,x)*(k*n*tt + p) = 0 - 5 - - 2 - 5*diff(tt,x)*k *n*tt -diff(p,x)*(k*n*tt - p) + diff(q,t) + diff(q,x)*u + ---------------------- - 2*m - - + diff(u,x)*q = 0 - - -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) - (u(j,m + 1) + u(j,m)))))/4 - - -1 - ht *(n(j,m + 1) - n(j,m) - n(j + 1,m) + n(j + 1,m + 1)) - + ---------------------------------------------------------- = 0 - 2 - - - -1 -( - hx *(n(j,m + 1)*tt(j,m + 1) + n(j,m)*tt(j,m) - n(j + 1,m)*tt(j + 1,m) - - -1 -1 - - n(j + 1,m + 1)*tt(j + 1,m + 1))*k)/2 + (hx *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)))*hx*m - - + 2*(p(j + 1,m + 1) + p(j + 1,m) - (p(j,m + 1) + p(j,m)))*ht + ( - - (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)))*ht*m) - - )/4 = 0 - - - -1 -(hx *((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)/4 - - -1 -1 -+ (hx *ht *(((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) - (q(j,m + 1) + q(j,m))))*ht + - - 3*((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)))*hx*k + 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)) - - )*ht*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 - - + 2*(5*p(j + 1,m + 1)*u(j + 1,m + 1) + 5*p(j + 1,m)*u(j + 1,m) - - - 5*p(j,m + 1)*u(j,m + 1) - 5*p(j,m)*u(j,m) + 2*q(j + 1,m + 1) - - + 2*q(j + 1,m) - 2*q(j,m + 1) - 2*q(j,m))))/20 - - -1 - ht *(p(j,m + 1) - p(j,m) - p(j + 1,m) + p(j + 1,m + 1)) - + ---------------------------------------------------------- = 0 - 2 - - - -1 -( - hx *(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)) - 2*( - - q(j + 1,m + 1)*u(j + 1,m + 1) + q(j + 1,m)*u(j + 1,m) - - - q(j,m + 1)*u(j,m + 1) - q(j,m)*u(j,m)))*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))*(tt(j + 1,m) - tt(j,m))) - - 2 - *k - 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))/(8*m) - - -1 - ht *(q(j,m + 1) - q(j,m) - q(j + 1,m) + q(j + 1,m + 1)) - + ---------------------------------------------------------- = 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; - -maxorder x=4; - - -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: - - 2 -hx - - 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] - - - [ hx 2*sin(ax)*ht*( - i*cos(ax) + sin(ax))] -h0 := [ ] - [ - 2*i*sin(ax)*ht hx ] - - - [ 2*sin(ax)*ht*( - i*cos(ax) + sin(ax)) ] - [ 1 ---------------------------------------] - [ hx ] -bb := [ ] - [ - 2*i*sin(ax)*ht ] - [------------------- 1 ] - [ hx ] - - -bb:=denotemat bb; - - - [ 1 ai12*i + ar12] -bb := [ ] - [ai21*i 1 ] - - -factor lam; - - -pol:=charpol bb; - - - 2 -pol := lam - 2*lam + ai12*ai21 - i*ai21*ar12 + 1 - -prdenot; - - - 2 - 2*sin(ax) *ht -ar12 := --------------- - hx - - - 2*cos(ax)*sin(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 := [cos(ax)*a*ht - i*sin(ax)*a*ht - a*ht + hx] - - - [ cos(ax)*a*ht - i*sin(ax)*a*ht - a*ht + hx ] -bb := [-------------------------------------------] - [ hx ] - - -bb:=denotemat bb; - - -bb := [ai11*i + 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 - - - [ i*sin(ax)*ht*r0 i*sin(ax)*ht*u0 + hx] -h1 := [ ] - [2*r0*(i*sin(ax)*ht*u0 + hx) 0 ] - - -h0 := mat((0,hx), - - 2 2 - (cos(ax)*c0 *ht - i*sin(ax)*c0 *ht + 2*hx*r0, - - 2 - c0 *ht*( - cos(ax) - i*sin(ax)))) - - - 2 2 - - i*cos(ax)*c0 *ht - sin(ax)*c0 *ht - 2*i*hx*r0 -bb := mat((--------------------------------------------------, - 2*r0*(sin(ax)*ht*u0 - i*hx) - - 2 - c0 *ht*(i*cos(ax) - sin(ax)) - ------------------------------), - 2*r0*(sin(ax)*ht*u0 - i*hx) - - 2 2 - sin(ax)*ht*(i*cos(ax)*c0 *ht + sin(ax)*c0 *ht + 2*i*hx*r0) - (------------------------------------------------------------,( - 2 2 2 2 - 2*(sin(ax) *ht *u0 - 2*i*sin(ax)*ht*hx*u0 - hx ) - - 2 2 2 2 2 - - i*cos(ax)*sin(ax)*c0 *ht + sin(ax) *c0 *ht - - 2 - - 2*i*sin(ax)*ht*hx*u0 - 2*hx )/(2 - - 2 2 2 2 - *(sin(ax) *ht *u0 - 2*i*sin(ax)*ht*hx*u0 - hx )))) - - -bb:=denotemat bb; - - - [ai11*i + ar11 ai12*i + ar12] -bb := [ ] - [ai21*i + ar21 ai22*i + ar22] - - -pol:=charpol bb; - - - 2 -pol := lam + lam*( - i*ai11 - i*ai22 - ar11 - ar22) - ai11*ai22 + i*ai11*ar22 - - + ai12*ai21 - i*ai12*ar21 - i*ai21*ar12 + i*ai22*ar11 + ar11*ar22 - - - ar12*ar21 - -prdenot; - - - 2 - - c0 -ar11 := --------- - 2*r0*u0 - - 2 2 - - cos(ax)*c0 *ht*u0 - c0 *hx - 2*hx*r0*u0 -ai11 := -------------------------------------------- - 2*r0*u0*(sin(ax)*ht*u0 - hx*i) - - 2 - - c0 -ar12 := --------- - 2*r0*u0 - - 2 - c0 *(cos(ax)*ht*u0 - hx) -ai12 := -------------------------------- - 2*r0*u0*(sin(ax)*ht*u0 - hx*i) - - 2 2 2 - sin(ax) *c0 *ht -ar21 := ---------------------------- - 2 2 2 2 - 2*(sin(ax) *ht *u0 - hx ) - - 2 2 3 2 2 2 -ai21 := (sin(ax)*ht*(cos(ax)*sin(ax) *c0 *ht *u0 - cos(ax)*c0 *ht*hx - - 2 2 2 2 2 2 3 - + 2*sin(ax) *c0 *ht *hx*u0 + 2*sin(ax) *ht *hx*r0*u0 - 2*hx *r0))/ - - 4 4 4 3 3 3 2 2 2 2 - (2*(sin(ax) *ht *u0 - 2*sin(ax) *ht *hx*i*u0 - 2*sin(ax) *ht *hx *u0 - - 3 4 - + 2*sin(ax)*ht*hx *i*u0 + hx )) - - 2 2 2 2 - sin(ax) *c0 *ht - 2*hx -ar22 := ---------------------------- - 2 2 2 2 - 2*(sin(ax) *ht *u0 - hx ) - - 2 2 3 2 2 2 -ai22 := (sin(ax)*ht*( - cos(ax)*sin(ax) *c0 *ht *u0 + cos(ax)*c0 *ht*hx - - 2 2 2 2 2 3 3 - + 2*sin(ax) *c0 *ht *hx*u0 - 2*sin(ax) *ht *hx*u0 - 2*hx *u0))/(2* - - 4 4 4 3 3 3 2 2 2 2 - (sin(ax) *ht *u0 - 2*sin(ax) *ht *hx*i*u0 - 2*sin(ax) *ht *hx *u0 - - 3 4 - + 2*sin(ax)*ht*hx *i*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 - - - [ hx - i*sin(ax)*c*ht ] -h1 := [ ] - [sin(ax)*c*ht*( - i*cos(ax) - sin(ax)) hx*(cos(ax) - i*sin(ax))] - - - [ hx i*sin(ax)*c*ht ] -h0 := [ ] - [sin(ax)*c*ht*(i*cos(ax) + sin(ax)) hx*(cos(ax) - i*sin(ax))] - - - [ 2 2 2 2 ] - [ - sin(ax) *c *ht + hx 2*i*sin(ax)*c*ht*hx ] - [-------------------------- ----------------------- ] - [ 2 2 2 2 2 2 2 2 ] - [ sin(ax) *c *ht + hx sin(ax) *c *ht + hx ] -bb := [ ] - [ 2 2 2 2 ] - [ 2*i*sin(ax)*c*ht*hx - 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 ai12*i] -bb := [ ] - [ai21*i 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*ht*hx -ai12 := ----------------------- - 2 2 2 2 - sin(ax) *c *ht + hx - - 2*sin(ax)*c*ht*hx -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 - - - cos(ax)*c*ht*hy + cos(ay)*c*ht*hx - c*ht*hx - c*ht*hy + hx*hy -bb := mat((---------------------------------------------------------------, - hx*hy - - - i*sin(ax)*c*ht - i*sin(ay)*c*ht - -------------------,-------------------), - hx hy - - - i*sin(ax)*c*ht cos(ax)*c*ht - c*ht + hx - (-------------------,--------------------------,0), - hx hx - - - i*sin(ay)*c*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*ht*hx - - + 2*c*ht*hx + 2*c*ht*hy - 3*hx*hy) + lam*( - - 2 2 2 2 - 3*cos(ax)*cos(ay)*c *ht *hx*hy - 3*cos(ax)*c *ht *hx*hy - - 2 2 2 2 2 2 2 2 - - 2*cos(ax)*c *ht *hy + 4*cos(ax)*c*ht*hx*hy + cos(ay) *c *ht *hx - - 2 2 2 2 2 - - 2*cos(ay)*c *ht *hx - 3*cos(ay)*c *ht *hx*hy - - 2 2 2 2 2 2 2 2 - + 4*cos(ay)*c*ht*hx *hy + sin(ay) *c *ht *hx + c *ht *hx - - 2 2 2 2 2 2 2 - + 3*c *ht *hx*hy + 2*c *ht *hy - 4*c*ht*hx *hy - 4*c*ht*hx*hy - - 2 2 2 3 3 - + 3*hx *hy ) - cos(ax)*cos(ay) *c *ht *hx - - 3 3 3 3 - + 2*cos(ax)*cos(ay)*c *ht *hx + 2*cos(ax)*cos(ay)*c *ht *hy - - 2 2 2 3 3 - - 3*cos(ax)*cos(ay)*c *ht *hx*hy - cos(ax)*sin(ay) *c *ht *hx - - 3 3 3 3 2 2 - - cos(ax)*c *ht *hx - 2*cos(ax)*c *ht *hy + 3*cos(ax)*c *ht *hx*hy - - 2 2 2 2 2 3 3 - + 2*cos(ax)*c *ht *hy - 2*cos(ax)*c*ht*hx*hy + cos(ay) *c *ht *hx - - 2 2 2 2 3 3 3 3 - - cos(ay) *c *ht *hx - 2*cos(ay)*c *ht *hx - 2*cos(ay)*c *ht *hy - - 2 2 2 2 2 2 - + 2*cos(ay)*c *ht *hx + 3*cos(ay)*c *ht *hx*hy - 2*cos(ay)*c*ht*hx *hy - - 2 3 3 2 2 2 2 3 3 3 3 - + sin(ay) *c *ht *hx - sin(ay) *c *ht *hx + c *ht *hx + 2*c *ht *hy - - 2 2 2 2 2 2 2 2 2 - - c *ht *hx - 3*c *ht *hx*hy - 2*c *ht *hy + 2*c*ht*hx *hy - - 2 2 2 2 2 - + 2*c*ht*hx*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 2 2 -pol := lam + lam *(4*s1 *cap1 + 4*s2 *cap2 - 3) + lam*(12*s1 *s2 *cap1*cap2 - - 2 2 2 2 2 2 - + 4*s1 *cap1 - 8*s1 *cap1 + 4*s2 *cap2 - 8*s2 *cap2 + 3) - - 2 2 2 2 2 2 2 2 - + 8*s1 *s2 *cap1 *cap2 + 8*s1 *s2 *cap1*cap2 - 12*s1 *s2 *cap1*cap2 - - 2 2 2 2 2 2 - - 4*s1 *cap1 + 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 2 2 -pol := lam + lam *(4*s1 *cap1 + 4*s2 *cap2 - 3) + lam*(12*s1 *s2 *cap1*cap2 - - 2 2 2 2 2 2 - + 4*s1 *cap1 - 8*s1 *cap1 + 4*s2 *cap2 - 8*s2 *cap2 + 3) - - 2 2 2 2 2 2 2 2 - + 8*s1 *s2 *cap1 *cap2 + 8*s1 *s2 *cap1*cap2 - 12*s1 *s2 *cap1*cap2 - - 2 2 2 2 2 2 - - 4*s1 *cap1 + 4*s1 *cap1 - 4*s2 *cap2 + 4*s2 *cap2 - 1 - -pol1:=hurw pol; - - - 3 2 2 2 2 2 2 -pol1 := 8*lam *s1 *s2 *cap1*cap2*(cap1 + cap2) + 8*lam *( - 3*s1 *s2 *cap1 *cap2 - - 2 2 2 2 2 2 2 2 2 - - 3*s1 *s2 *cap1*cap2 + 3*s1 *s2 *cap1*cap2 + s1 *cap1 + s2 *cap2 - - 2 2 2 2 2 2 - ) + 8*lam*(3*s1 *s2 *cap1 *cap2 + 3*s1 *s2 *cap1*cap2 - - 2 2 2 2 2 2 2 - - 6*s1 *s2 *cap1*cap2 - 2*s1 *cap1 + 2*s1 *cap1 - 2*s2 *cap2 - - 2 2 2 2 2 2 2 - + 2*s2 *cap2) + 8*( - s1 *s2 *cap1 *cap2 - s1 *s2 *cap1*cap2 - - 2 2 2 2 2 2 2 - + 3*s1 *s2 *cap1*cap2 + s1 *cap1 - 2*s1 *cap1 + s2 *cap2 - - 2 - - 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 2 2 -cpr00 := 8*s1 *s2 *cap1 *cap2 + 8*s1 *s2 *cap1*cap2 - 12*s1 *s2 *cap1*cap2 - - 2 2 2 2 2 2 - - 4*s1 *cap1 + 4*s1 *cap1 - 4*s2 *cap2 + 4*s2 *cap2 - 1 - -cpr01 := - - 2 2 2 2 2 2 2 2 -12*s1 *s2 *cap1*cap2 + 4*s1 *cap1 - 8*s1 *cap1 + 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 - - -bb := mat((( - i*cos(ax)*sin(ax)*ht*hx*w - i*cos(ax)*sin(ay)*ht*hx*v - - - i*cos(ay)*sin(ax)*ht*hx*w - i*cos(ay)*sin(ay)*ht*hx*v - - 2 2 2 2 2 2 2 - - 2*sin(ax) *ht *w - 4*sin(ax)*sin(ay)*ht *v*w - 2*sin(ay) *ht *v - - 2 2 - + hx )/hx ,(sin(ax)*ht*p*( - i*cos(ax)*hx - i*cos(ay)*hx - - 2 - - 4*sin(ax)*ht*w - 4*sin(ay)*ht*v))/hx ,(ht*p*( - - - i*cos(ax)*sin(ay)*hx - 4*i*cos(ay)*sin(ay)*ht*v - - 2 - - i*cos(ay)*sin(ay)*hx - 4*sin(ax)*sin(ay)*ht*w - 2*ht*v))/hx - - 2 2 2 - - 2*ht *(sin(ax) + sin(ay) ) - ,--------------------------------), - 2 - hx - - (0,( - i*cos(ax)*sin(ax)*ht*hx*w - i*cos(ax)*sin(ay)*ht*hx*v - - - i*cos(ay)*sin(ax)*ht*hx*w - i*cos(ay)*sin(ay)*ht*hx*v - - 2 2 2 2 2 2 - - 2*sin(ax) *c *ht - 2*sin(ax) *ht *w - - 2 2 2 2 2 2 - - 4*sin(ax)*sin(ay)*ht *v*w - 2*sin(ay) *ht *v + hx )/hx , - - 2 2 - - 2*sin(ax)*sin(ay)*c *ht - -----------------------------,(sin(ax)*ht*( - i*cos(ax)*hx - 2 - hx - - 2 - - i*cos(ay)*hx - 4*sin(ax)*ht*w - 4*sin(ay)*ht*v))/(hx *p)), - - 2 2 - - 2*sin(ax)*sin(ay)*c *ht - (0,-----------------------------,( - i*cos(ax)*sin(ax)*ht*hx*w - 2 - hx - - - i*cos(ax)*sin(ay)*ht*hx*v - i*cos(ay)*sin(ax)*ht*hx*w - - 2 2 2 - - i*cos(ay)*sin(ay)*ht*hx*v - 2*sin(ax) *ht *w - - 2 2 2 2 - - 4*sin(ax)*sin(ay)*ht *v*w - 2*sin(ay) *c *ht - - 2 2 2 2 2 - - 2*sin(ay) *ht *v + hx )/hx ,(ht*( - 2*cos(ax)*cos(ay)*ht*w - - - 2*i*cos(ax)*sin(ay)*ht*w - i*cos(ax)*sin(ay)*hx - - - 2*i*cos(ay)*sin(ax)*ht*w - i*cos(ay)*sin(ay)*hx - - 2 2 - - 2*sin(ax)*sin(ay)*ht*w - 4*sin(ay) *ht*v))/(hx *p)), - - 2 - (0,(sin(ax)*c *ht*p*( - i*cos(ax)*hx - i*cos(ay)*hx - 4*sin(ax)*ht*w - - 2 2 - - 4*sin(ay)*ht*v))/hx ,(sin(ay)*c *ht*p*( - i*cos(ax)*hx - - 2 - - i*cos(ay)*hx - 4*sin(ax)*ht*w - 4*sin(ay)*ht*v))/hx ,( - - - i*cos(ax)*sin(ax)*ht*hx*w - i*cos(ax)*sin(ay)*ht*hx*v - - - i*cos(ay)*sin(ax)*ht*hx*w - i*cos(ay)*sin(ay)*ht*hx*v - - 2 2 2 2 2 2 - - 2*sin(ax) *c *ht - 2*sin(ax) *ht *w - - 2 2 2 2 - - 4*sin(ax)*sin(ay)*ht *v*w - 2*sin(ay) *c *ht - - 2 2 2 2 2 - - 2*sin(ay) *ht *v + 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; - - - [ai11*i + ar11 ai12*i + ar12 ai13*i + ar13 ar14 ] - [ ] - [ 0 ai22*i + ar22 ar23 ai24*i + ar24] -bb := [ ] - [ 0 ar32 ai33*i + ar33 ai34*i + ar34] - [ ] - [ 0 ai42*i + ar42 ai43*i + ar43 ai44*i + 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 - ar33 - ar44) + - - 2 - lam *( - ai11*ai22 - ai11*ai33 - ai11*ai44 + i*ai11*ar22 + i*ai11*ar33 - - + i*ai11*ar44 - ai22*ai33 - ai22*ai44 + i*ai22*ar11 + i*ai22*ar33 - - + i*ai22*ar44 + ai24*ai42 - i*ai24*ar42 - ai33*ai44 + i*ai33*ar11 - - + i*ai33*ar22 + i*ai33*ar44 + ai34*ai43 - i*ai34*ar43 - - - i*ai42*ar24 - i*ai43*ar34 + i*ai44*ar11 + i*ai44*ar22 - - + i*ai44*ar33 + ar11*ar22 + ar11*ar33 + ar11*ar44 + ar22*ar33 - - + ar22*ar44 - ar23*ar32 - ar24*ar42 + ar33*ar44 - ar34*ar43) + lam - - *(i*ai11*ai22*ai33 + i*ai11*ai22*ai44 + ai11*ai22*ar33 + ai11*ai22*ar44 - - - i*ai11*ai24*ai42 - ai11*ai24*ar42 + i*ai11*ai33*ai44 - - + ai11*ai33*ar22 + ai11*ai33*ar44 - i*ai11*ai34*ai43 - ai11*ai34*ar43 - - - ai11*ai42*ar24 - ai11*ai43*ar34 + ai11*ai44*ar22 + ai11*ai44*ar33 - - - i*ai11*ar22*ar33 - i*ai11*ar22*ar44 + i*ai11*ar23*ar32 - - + i*ai11*ar24*ar42 - i*ai11*ar33*ar44 + i*ai11*ar34*ar43 - - + i*ai22*ai33*ai44 + ai22*ai33*ar11 + ai22*ai33*ar44 - - - i*ai22*ai34*ai43 - ai22*ai34*ar43 - ai22*ai43*ar34 + ai22*ai44*ar11 - - + ai22*ai44*ar33 - i*ai22*ar11*ar33 - i*ai22*ar11*ar44 - - - i*ai22*ar33*ar44 + i*ai22*ar34*ar43 - i*ai24*ai33*ai42 - - - ai24*ai33*ar42 - ai24*ai42*ar11 - ai24*ai42*ar33 + ai24*ai43*ar32 - - + i*ai24*ar11*ar42 - i*ai24*ar32*ar43 + i*ai24*ar33*ar42 - - - ai33*ai42*ar24 + ai33*ai44*ar11 + ai33*ai44*ar22 - i*ai33*ar11*ar22 - - - i*ai33*ar11*ar44 - i*ai33*ar22*ar44 + i*ai33*ar24*ar42 - - + ai34*ai42*ar23 - ai34*ai43*ar11 - ai34*ai43*ar22 + i*ai34*ar11*ar43 - - + i*ai34*ar22*ar43 - i*ai34*ar23*ar42 + i*ai42*ar11*ar24 - - - i*ai42*ar23*ar34 + i*ai42*ar24*ar33 + i*ai43*ar11*ar34 - - + i*ai43*ar22*ar34 - i*ai43*ar24*ar32 - i*ai44*ar11*ar22 - - - i*ai44*ar11*ar33 - i*ai44*ar22*ar33 + i*ai44*ar23*ar32 - - - ar11*ar22*ar33 - ar11*ar22*ar44 + ar11*ar23*ar32 + ar11*ar24*ar42 - - - ar11*ar33*ar44 + ar11*ar34*ar43 - ar22*ar33*ar44 + ar22*ar34*ar43 - - + ar23*ar32*ar44 - ar23*ar34*ar42 - ar24*ar32*ar43 + ar24*ar33*ar42) - - + ai11*ai22*ai33*ai44 - i*ai11*ai22*ai33*ar44 - ai11*ai22*ai34*ai43 - - + i*ai11*ai22*ai34*ar43 + i*ai11*ai22*ai43*ar34 - i*ai11*ai22*ai44*ar33 - - - ai11*ai22*ar33*ar44 + ai11*ai22*ar34*ar43 - ai11*ai24*ai33*ai42 - - + i*ai11*ai24*ai33*ar42 + i*ai11*ai24*ai42*ar33 - i*ai11*ai24*ai43*ar32 - - - ai11*ai24*ar32*ar43 + ai11*ai24*ar33*ar42 + i*ai11*ai33*ai42*ar24 - - - i*ai11*ai33*ai44*ar22 - ai11*ai33*ar22*ar44 + ai11*ai33*ar24*ar42 - - - i*ai11*ai34*ai42*ar23 + i*ai11*ai34*ai43*ar22 + ai11*ai34*ar22*ar43 - - - ai11*ai34*ar23*ar42 - ai11*ai42*ar23*ar34 + ai11*ai42*ar24*ar33 - - + ai11*ai43*ar22*ar34 - ai11*ai43*ar24*ar32 - ai11*ai44*ar22*ar33 - - + ai11*ai44*ar23*ar32 + i*ai11*ar22*ar33*ar44 - i*ai11*ar22*ar34*ar43 - - - i*ai11*ar23*ar32*ar44 + i*ai11*ar23*ar34*ar42 + i*ai11*ar24*ar32*ar43 - - - i*ai11*ar24*ar33*ar42 - i*ai22*ai33*ai44*ar11 - ai22*ai33*ar11*ar44 - - + i*ai22*ai34*ai43*ar11 + ai22*ai34*ar11*ar43 + ai22*ai43*ar11*ar34 - - - ai22*ai44*ar11*ar33 + i*ai22*ar11*ar33*ar44 - i*ai22*ar11*ar34*ar43 - - + i*ai24*ai33*ai42*ar11 + ai24*ai33*ar11*ar42 + ai24*ai42*ar11*ar33 - - - ai24*ai43*ar11*ar32 + i*ai24*ar11*ar32*ar43 - i*ai24*ar11*ar33*ar42 - - + ai33*ai42*ar11*ar24 - ai33*ai44*ar11*ar22 + i*ai33*ar11*ar22*ar44 - - - i*ai33*ar11*ar24*ar42 - ai34*ai42*ar11*ar23 + ai34*ai43*ar11*ar22 - - - i*ai34*ar11*ar22*ar43 + i*ai34*ar11*ar23*ar42 + i*ai42*ar11*ar23*ar34 - - - i*ai42*ar11*ar24*ar33 - i*ai43*ar11*ar22*ar34 + i*ai43*ar11*ar24*ar32 - - + i*ai44*ar11*ar22*ar33 - i*ai44*ar11*ar23*ar32 + ar11*ar22*ar33*ar44 - - - ar11*ar22*ar34*ar43 - ar11*ar23*ar32*ar44 + ar11*ar23*ar34*ar42 - - + ar11*ar24*ar32*ar43 - ar11*ar24*ar33*ar42 - -denotid cp; - - -pol:=denotepol pol; - - - 4 3 2 -pol := lam + lam *(cpi03*i + cpr03) + lam *(cpi02*i + cpr02) - - + lam*(cpi01*i + cpr01) + cpi00*i + 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 *(cpi03 + 2*cpr02 + cpr03 ) - - 5 - + 2*lam *(cpi02*cpi03 + cpr01 + cpr02*cpr03) - - 4 2 2 - + lam *(2*cpi01*cpi03 + cpi02 + 2*cpr00 + 2*cpr01*cpr03 + cpr02 ) - - 3 - + 2*lam *(cpi00*cpi03 + cpi01*cpi02 + cpr00*cpr03 + cpr01*cpr02) - - 2 2 2 - + lam *(2*cpi00*cpi02 + cpi01 + 2*cpr00*cpr02 + cpr01 ) - - 2 2 - + 2*lam*(cpi00*cpi01 + cpr00*cpr01) + cpi00 + cpr00 - -denotid rp; - - -pol:=denotepol pol; - - - 8 7 6 5 4 3 -pol := lam + lam *rpr07 + lam *rpr06 + lam *rpr05 + lam *rpr04 + lam *rpr03 - - 2 - + 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*k3 *p*(s1*k1 + s2*k2) -ar42 := ------------------------------- - r1 - - 2 - - s1*k3 *p*(c1 + c2) -ai42 := ----------------------- - r1 - - 2 - - 4*s2*k3 *p*(s1*k1 + s2*k2) -ar43 := ------------------------------- - r1 - - 2 - - s2*k3 *p*(c1 + c2) -ai43 := ----------------------- - r1 - - 2 2 2 2 2 2 2 2 -ar44 := - 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 := ai11*ai22*ai33*ai44 - ai11*ai22*ai34*ai43 - ai11*ai22*ar33*ar44 - - + ai11*ai22*ar34*ar43 - ai11*ai24*ai33*ai42 - ai11*ai24*ar32*ar43 - - + ai11*ai24*ar33*ar42 - ai11*ai33*ar22*ar44 + ai11*ai33*ar24*ar42 - - + ai11*ai34*ar22*ar43 - ai11*ai34*ar23*ar42 - ai11*ai42*ar23*ar34 - - + ai11*ai42*ar24*ar33 + ai11*ai43*ar22*ar34 - ai11*ai43*ar24*ar32 - - - ai11*ai44*ar22*ar33 + ai11*ai44*ar23*ar32 - ai22*ai33*ar11*ar44 - - + ai22*ai34*ar11*ar43 + ai22*ai43*ar11*ar34 - ai22*ai44*ar11*ar33 - - + ai24*ai33*ar11*ar42 + ai24*ai42*ar11*ar33 - ai24*ai43*ar11*ar32 - - + ai33*ai42*ar11*ar24 - ai33*ai44*ar11*ar22 - ai34*ai42*ar11*ar23 - - + ai34*ai43*ar11*ar22 + ar11*ar22*ar33*ar44 - ar11*ar22*ar34*ar43 - - - ar11*ar23*ar32*ar44 + ar11*ar23*ar34*ar42 + ar11*ar24*ar32*ar43 - - - ar11*ar24*ar33*ar42 - -cpi00 := - ai11*ai22*ai33*ar44 + ai11*ai22*ai34*ar43 + ai11*ai22*ai43*ar34 - - - ai11*ai22*ai44*ar33 + ai11*ai24*ai33*ar42 + ai11*ai24*ai42*ar33 - - - ai11*ai24*ai43*ar32 + ai11*ai33*ai42*ar24 - ai11*ai33*ai44*ar22 - - - ai11*ai34*ai42*ar23 + ai11*ai34*ai43*ar22 + ai11*ar22*ar33*ar44 - - - ai11*ar22*ar34*ar43 - ai11*ar23*ar32*ar44 + ai11*ar23*ar34*ar42 - - + ai11*ar24*ar32*ar43 - ai11*ar24*ar33*ar42 - ai22*ai33*ai44*ar11 - - + ai22*ai34*ai43*ar11 + ai22*ar11*ar33*ar44 - ai22*ar11*ar34*ar43 - - + ai24*ai33*ai42*ar11 + ai24*ar11*ar32*ar43 - ai24*ar11*ar33*ar42 - - + ai33*ar11*ar22*ar44 - ai33*ar11*ar24*ar42 - ai34*ar11*ar22*ar43 - - + ai34*ar11*ar23*ar42 + ai42*ar11*ar23*ar34 - ai42*ar11*ar24*ar33 - - - ai43*ar11*ar22*ar34 + ai43*ar11*ar24*ar32 + ai44*ar11*ar22*ar33 - - - ai44*ar11*ar23*ar32 - -cpr01 := ai11*ai22*ar33 + ai11*ai22*ar44 - ai11*ai24*ar42 + ai11*ai33*ar22 - - + ai11*ai33*ar44 - ai11*ai34*ar43 - ai11*ai42*ar24 - ai11*ai43*ar34 - - + ai11*ai44*ar22 + ai11*ai44*ar33 + ai22*ai33*ar11 + ai22*ai33*ar44 - - - ai22*ai34*ar43 - ai22*ai43*ar34 + ai22*ai44*ar11 + ai22*ai44*ar33 - - - ai24*ai33*ar42 - ai24*ai42*ar11 - ai24*ai42*ar33 + ai24*ai43*ar32 - - - ai33*ai42*ar24 + ai33*ai44*ar11 + ai33*ai44*ar22 + ai34*ai42*ar23 - - - ai34*ai43*ar11 - ai34*ai43*ar22 - ar11*ar22*ar33 - ar11*ar22*ar44 - - + ar11*ar23*ar32 + ar11*ar24*ar42 - ar11*ar33*ar44 + ar11*ar34*ar43 - - - ar22*ar33*ar44 + ar22*ar34*ar43 + ar23*ar32*ar44 - ar23*ar34*ar42 - - - ar24*ar32*ar43 + ar24*ar33*ar42 - -cpi01 := ai11*ai22*ai33 + ai11*ai22*ai44 - ai11*ai24*ai42 + ai11*ai33*ai44 - - - ai11*ai34*ai43 - ai11*ar22*ar33 - ai11*ar22*ar44 + ai11*ar23*ar32 - - + ai11*ar24*ar42 - ai11*ar33*ar44 + ai11*ar34*ar43 + ai22*ai33*ai44 - - - ai22*ai34*ai43 - ai22*ar11*ar33 - ai22*ar11*ar44 - ai22*ar33*ar44 - - + ai22*ar34*ar43 - ai24*ai33*ai42 + ai24*ar11*ar42 - ai24*ar32*ar43 - - + ai24*ar33*ar42 - ai33*ar11*ar22 - ai33*ar11*ar44 - ai33*ar22*ar44 - - + ai33*ar24*ar42 + ai34*ar11*ar43 + ai34*ar22*ar43 - ai34*ar23*ar42 - - + ai42*ar11*ar24 - ai42*ar23*ar34 + ai42*ar24*ar33 + ai43*ar11*ar34 - - + ai43*ar22*ar34 - ai43*ar24*ar32 - ai44*ar11*ar22 - ai44*ar11*ar33 - - - ai44*ar22*ar33 + ai44*ar23*ar32 - -cpr02 := - ai11*ai22 - ai11*ai33 - ai11*ai44 - ai22*ai33 - ai22*ai44 - - + ai24*ai42 - ai33*ai44 + ai34*ai43 + ar11*ar22 + ar11*ar33 - - + ar11*ar44 + ar22*ar33 + ar22*ar44 - ar23*ar32 - ar24*ar42 - - + ar33*ar44 - ar34*ar43 - -cpi02 := ai11*ar22 + ai11*ar33 + ai11*ar44 + ai22*ar11 + ai22*ar33 + ai22*ar44 - - - ai24*ar42 + ai33*ar11 + ai33*ar22 + ai33*ar44 - ai34*ar43 - - - ai42*ar24 - ai43*ar34 + ai44*ar11 + ai44*ar22 + ai44*ar33 - -cpr03 := - (ar11 + ar22 + ar33 + ar44) - -cpi03 := - (ai11 + ai22 + ai33 + ai44) - - 2 2 -rpr00 := cpi00 + cpr00 - -rpr01 := 2*(cpi00*cpi01 + cpr00*cpr01) - - 2 2 -rpr02 := 2*cpi00*cpi02 + cpi01 + 2*cpr00*cpr02 + cpr01 - -rpr03 := 2*(cpi00*cpi03 + cpi01*cpi02 + cpr00*cpr03 + cpr01*cpr02) - - 2 2 -rpr04 := 2*cpi01*cpi03 + cpi02 + 2*cpr00 + 2*cpr01*cpr03 + cpr02 - -rpr05 := 2*(cpi02*cpi03 + cpr01 + cpr02*cpr03) - - 2 2 -rpr06 := cpi03 + 2*cpr02 + cpr03 - -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 - (ai*i + ar) - -x2:=lam-(br+i*bi); - - -x2 := lam - (bi*i + br) - -x3:=lam-(cr+i*ci); - - -x3 := lam - (ci*i + cr) - -hurwitzp x1; - -Necessary and sufficient conditions are: - - ar > 0 - - -cond - - -% Example H.2 - -x:=hurw(x0*x1); - - -x := 2*lam*( - ai*i - ar + 1) + 2*(ai*i + ar + 1) - -hurwitzp x; - -Necessary and sufficient conditions are: - 2 2 -4*( - ai - ar + 1) > 0 - - -cond - - -% Example H.3 - -x:=(x1*x2); - - - 2 -x := lam - lam*(ai*i + ar + bi*i + br) - ai*bi + ai*br*i + ar*bi*i + ar*br - -hurwitzp x; - -Necessary and sufficient conditions are: - - (ar + br) > 0 - - 2 2 2 2 -ar*br*(ai - 2*ai*bi + ar + 2*ar*br + bi + br ) > 0 - - -cond - - -% Example H.4 - -x:=(x1*x2*x3); - - - 3 2 -x := lam - lam *(ai*i + ar + bi*i + br + ci*i + cr) + lam*( - ai*bi + ai*br*i - - - ai*ci + ai*cr*i + ar*bi*i + ar*br + ar*ci*i + ar*cr - bi*ci + bi*cr*i - - + br*ci*i + br*cr) + ai*bi*ci*i + ai*bi*cr + ai*br*ci - ai*br*cr*i - - + ar*bi*ci - ar*bi*cr*i - ar*br*ci*i - ar*br*cr - -hurwitzp x; - -Necessary and sufficient conditions are: - - (ar + br + cr) > 0 - - 2 2 3 3 -ai *ar*br + ai *ar*cr - 2*ai*ar*bi*br - 2*ai*ar*ci*cr + ar *br + ar *cr - - 2 2 2 2 2 2 3 2 - + 2*ar *br + 4*ar *br*cr + 2*ar *cr + ar*bi *br + ar*br + 4*ar*br *cr - - 2 2 3 2 3 - + 4*ar*br*cr + ar*ci *cr + ar*cr + bi *br*cr - 2*bi*br*ci*cr + br *cr - - 2 2 2 3 - + 2*br *cr + br*ci *cr + br*cr > 0 - - 4 2 4 4 2 4 4 2 4 2 -ar*br*cr*( - ai *bi + 2*ai *bi*ci - ai *br - 2*ai *br*cr - ai *ci - ai *cr - - 3 3 3 2 3 2 3 - + 2*ai *bi - 2*ai *bi *ci + 2*ai *bi*br + 4*ai *bi*br*cr - - 3 2 3 2 3 2 3 - - 2*ai *bi*ci + 2*ai *bi*cr + 2*ai *br *ci + 4*ai *br*ci*cr - - 3 3 3 2 2 2 2 2 2 - + 2*ai *ci + 2*ai *ci*cr - 2*ai *ar *bi + 4*ai *ar *bi*ci - - 2 2 2 2 2 2 2 2 2 2 2 - - 2*ai *ar *br - 4*ai *ar *br*cr - 2*ai *ar *ci - 2*ai *ar *cr - - 2 2 2 2 2 - - 2*ai *ar*bi *br - 2*ai *ar*bi *cr + 4*ai *ar*bi*br*ci - - 2 2 3 2 2 - + 4*ai *ar*bi*ci*cr - 2*ai *ar*br - 6*ai *ar*br *cr - - 2 2 2 2 2 2 2 3 - - 2*ai *ar*br*ci - 6*ai *ar*br*cr - 2*ai *ar*ci *cr - 2*ai *ar*cr - - 2 4 2 3 2 2 2 2 2 - - ai *bi - 2*ai *bi *ci - 2*ai *bi *br - 2*ai *bi *br*cr - - 2 2 2 2 2 2 2 2 2 - + 6*ai *bi *ci - 2*ai *bi *cr - 2*ai *bi*br *ci - 8*ai *bi*br*ci*cr - - 2 3 2 2 2 4 2 3 - - 2*ai *bi*ci - 2*ai *bi*ci*cr - ai *br - 2*ai *br *cr - - 2 2 2 2 2 2 2 2 2 3 - - 2*ai *br *ci - 2*ai *br *cr - 2*ai *br*ci *cr - 2*ai *br*cr - - 2 4 2 2 2 2 4 2 3 2 2 - - ai *ci - 2*ai *ci *cr - ai *cr + 2*ai*ar *bi - 2*ai*ar *bi *ci - - 2 2 2 2 2 - + 2*ai*ar *bi*br + 4*ai*ar *bi*br*cr - 2*ai*ar *bi*ci - - 2 2 2 2 2 - + 2*ai*ar *bi*cr + 2*ai*ar *br *ci + 4*ai*ar *br*ci*cr - - 2 3 2 2 3 2 - + 2*ai*ar *ci + 2*ai*ar *ci*cr + 4*ai*ar*bi *cr + 4*ai*ar*bi *br*ci - - 2 2 2 - - 8*ai*ar*bi *ci*cr + 4*ai*ar*bi*br *cr - 8*ai*ar*bi*br*ci - - 2 2 3 - + 8*ai*ar*bi*br*cr + 4*ai*ar*bi*ci *cr + 4*ai*ar*bi*cr - - 3 2 3 - + 4*ai*ar*br *ci + 8*ai*ar*br *ci*cr + 4*ai*ar*br*ci - - 2 4 3 2 3 2 - + 4*ai*ar*br*ci*cr + 2*ai*bi *ci - 2*ai*bi *ci + 2*ai*bi *cr - - 2 2 2 2 3 - + 4*ai*bi *br *ci + 4*ai*bi *br*ci*cr - 2*ai*bi *ci - - 2 2 2 2 2 2 - - 2*ai*bi *ci*cr - 2*ai*bi*br *ci + 2*ai*bi*br *cr - - 2 3 4 2 2 - + 4*ai*bi*br*ci *cr + 4*ai*bi*br*cr + 2*ai*bi*ci + 4*ai*bi*ci *cr - - 4 4 3 2 3 - + 2*ai*bi*cr + 2*ai*br *ci + 4*ai*br *ci*cr + 2*ai*br *ci - - 2 2 4 2 4 4 2 4 - + 2*ai*br *ci*cr - ar *bi + 2*ar *bi*ci - ar *br - 2*ar *br*cr - - 4 2 4 2 3 2 3 2 3 - - ar *ci - ar *cr - 2*ar *bi *br - 2*ar *bi *cr + 4*ar *bi*br*ci - - 3 3 3 3 2 3 2 - + 4*ar *bi*ci*cr - 2*ar *br - 6*ar *br *cr - 2*ar *br*ci - - 3 2 3 2 3 3 2 4 2 3 - - 6*ar *br*cr - 2*ar *ci *cr - 2*ar *cr - ar *bi + 2*ar *bi *ci - - 2 2 2 2 2 2 2 2 2 2 2 - - 2*ar *bi *br - 6*ar *bi *br*cr - 2*ar *bi *ci - 2*ar *bi *cr - - 2 2 2 2 3 - + 2*ar *bi*br *ci + 8*ar *bi*br*ci*cr + 2*ar *bi*ci - - 2 2 2 4 2 3 2 2 2 - + 2*ar *bi*ci*cr - ar *br - 6*ar *br *cr - 2*ar *br *ci - - 2 2 2 2 2 2 3 2 4 - - 10*ar *br *cr - 6*ar *br*ci *cr - 6*ar *br*cr - ar *ci - - 2 2 2 2 4 4 3 - - 2*ar *ci *cr - ar *cr - 2*ar*bi *cr + 4*ar*bi *ci*cr - - 2 2 2 2 2 2 - - 4*ar*bi *br *cr - 2*ar*bi *br*ci - 6*ar*bi *br*cr - - 2 2 2 3 2 3 - - 2*ar*bi *ci *cr - 2*ar*bi *cr + 4*ar*bi*br *ci*cr + 4*ar*bi*br*ci - - 2 4 3 2 3 2 - + 4*ar*bi*br*ci*cr - 2*ar*br *cr - 2*ar*br *ci - 6*ar*br *cr - - 2 2 2 3 4 2 2 - - 6*ar*br *ci *cr - 6*ar*br *cr - 2*ar*br*ci - 4*ar*br*ci *cr - - 4 4 2 4 2 3 3 3 2 - - 2*ar*br*cr - bi *ci - bi *cr + 2*bi *ci + 2*bi *ci*cr - - 2 2 2 2 2 2 2 2 2 3 - - 2*bi *br *ci - 2*bi *br *cr - 2*bi *br*ci *cr - 2*bi *br*cr - - 2 4 2 2 2 2 4 2 3 2 2 - - bi *ci - 2*bi *ci *cr - bi *cr + 2*bi*br *ci + 2*bi*br *ci*cr - - 4 2 4 2 3 2 3 3 2 4 - - br *ci - br *cr - 2*br *ci *cr - 2*br *cr - br *ci - - 2 2 2 2 4 - - 2*br *ci *cr - br *cr ) > 0 - - -cond - -clear x,x0,x1,x2,x3; - - - -%*********************************************************************** -%***** ***** -%***** T e s t Examples --- Module L I N B A N D ***** -%***** ***** -%*********************************************************************** - -on evallhseqp; - - % So both sides of equations evaluate. - -% 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=0.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.0e+000 - 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.0e+000 - 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.0 - arer=0.0 - 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=0.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.0e+000 - 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.0e+000 - 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.0 - arer=0.0 - 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=0.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.0e+000 - acof(n,2)=1.0 - acof(n,3)=0.0e+000 - 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.0e+000 - acof(n,2)=1.0 - acof(n,3)=0.0e+000 - 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.0 - arer=0.0 - 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=0.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.0e+000 - 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.0e+000 - 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.0 - arer=0.0 - 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 v(x..one),w(x..half); - - -iim aa, v, diff(v,t)=c*diff(w,x), - w, diff(w,t)=c*diff(v,x); - - -***************************** -***** Program ***** IIMET Ver 1.1.2 -***************************** - - 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 one grid point - Equation for w variable is integrated in half grid point - - Equations after Discretization Using IIM : - ========================================== - - 2*j - 1 2*j + 1 2*j + 1 2*j - 1 -((w(n,---------) - w(n,---------) - w(n + 1,---------) + w(n + 1,---------))*c - 2 2 2 2 - - *ht + 2*(v(n + 1,j) - v(n,j))*hx)/(2*ht*hx) = 0 - - -( - (v(n,j + 1) - v(n,j) - v(n + 1,j) + v(n + 1,j + 1))*c*ht - - 2*j + 1 2*j + 1 - + 2*(w(n + 1,---------) - w(n,---------))*hx)/(2*ht*hx) = 0 - 2 2 - - - -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 - - 2 -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 - - 2 -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 -a(1,1) := (2*v(n + 1,j)*hx - 2*v(n,j)*hx + w(n + 1,---------)*c*ht - 2 - - 2*j + 1 2*j - 1 - - w(n + 1,---------)*c*ht + w(n,---------)*c*ht - 2 2 - - 2*j + 1 - - w(n,---------)*c*ht)/(2*ht*hx) - 2 - -a(1,2):=aa(1,0); - - -a(1,2) := ( - v(n + 1,j + 1)*c*ht + v(n + 1,j)*c*ht - v(n,j + 1)*c*ht - - 2*j + 1 2*j + 1 - + v(n,j)*c*ht + 2*w(n + 1,---------)*hx - 2*w(n,---------)*hx)/(2*ht - 2 2 - - *hx) - -off prfourmat; - - -b:=ampmat a; - - - kx*hx -ax := ------- - 2 - - - [ 2 2 2 2 ] - [ - sin(ax) *c *ht + hx 2*i*sin(ax)*c*ht*hx ] - [-------------------------- ----------------------- ] - [ 2 2 2 2 2 2 2 2 ] - [ sin(ax) *c *ht + hx sin(ax) *c *ht + hx ] -b := [ ] - [ 2 2 2 2 ] - [ 2*i*sin(ax)*c*ht*hx - 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 *ht *hx + 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 2 2 2 - + 2*sin(ax) *c *ht *hx + hx )/(sin(ax) *c *ht + 2*sin(ax) *c *ht *hx - - 4 - + 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 *ht *hx + 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 2 2 2 - + 2*sin(ax) *c *ht *hx + hx )/(sin(ax) *c *ht + 2*sin(ax) *c *ht *hx - - 4 - + hx ) - -pol:=hurw num pol; - - -pol := - - 2 2 2 2 2 2 2 2 2 2 2 2 2 -4*lam *sin(ax) *c *ht *(sin(ax) *c *ht + hx ) + 4*hx *(sin(ax) *c *ht + hx ) - -hurwitzp pol; - -Necessary and sufficient conditions are: - 2 - hx ------------------ > 0 - 2 2 2 - sin(ax) *c *ht - - -cond - -bt:=tcon b; - - - [ 2 2 2 2 ] - [ - sin(ax) *c *ht + hx - 2*sin(ax)*c*ht*hx*i ] - [-------------------------- ------------------------ ] - [ 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*ht*hx*i - 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*i*sin(ax)*c*ht ] - [ 1 ------------------ ] - [ hx ] - [ ] -b := [ 2 2 2 2 ] - [ 2*i*sin(ax)*c*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 - - -cond - -bt:=tcon b; - - - [ - 2*sin(ax)*c*ht*i ] - [ 1 --------------------- ] - [ hx ] - [ ] -bt := [ 2 2 2 2 ] - [ - 2*sin(ax)*c*ht*i - 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 *ht *i ] -[------------------------- --------------------- ] -[ 2 3 ] -[ hx hx ] -[ ] -[ 3 3 3 4 4 4 2 2 2 2 4 ] -[ - 8*sin(ax) *c *ht *i 16*sin(ax) *c *ht - 4*sin(ax) *c *ht *hx + hx ] -[------------------------ --------------------------------------------------] -[ 3 4 ] -[ hx hx ] - - -bt*b-b*bt; - - -[ 3 3 3 ] -[ 16*sin(ax) *c *ht *i ] -[ 0 ----------------------] -[ 3 ] -[ hx ] -[ ] -[ 3 3 3 ] -[ - 16*sin(ax) *c *ht *i ] -[------------------------- 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 - i*sin(ax)*c*ht -b := mat((-------------------------------------------,-------------------, - hx hx - - - i*sin(ay)*c*ht - -------------------), - hx - - - i*sin(ax)*c*ht cos(ax)*c*ht - c*ht + hx - (-------------------,--------------------------,0), - hx hx - - - i*sin(ay)*c*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) - - 2 2 2 2 - + lam*hx*(3*cos(ax)*cos(ay)*c *ht - 5*cos(ax)*c *ht - - 2 2 - + 4*cos(ax)*c*ht*hx - 5*cos(ay)*c *ht + 4*cos(ay)*c*ht*hx - - 2 2 2 3 3 - + 7*c *ht - 8*c*ht*hx + 3*hx ) + 4*cos(ax)*cos(ay)*c *ht - - 2 2 3 3 2 2 - - 3*cos(ax)*cos(ay)*c *ht *hx - 4*cos(ax)*c *ht + 5*cos(ax)*c *ht *hx - - 2 3 3 2 2 - - 2*cos(ax)*c*ht*hx - 4*cos(ay)*c *ht + 5*cos(ay)*c *ht *hx - - 2 3 3 2 2 2 3 3 - - 2*cos(ay)*c*ht*hx + 4*c *ht - 7*c *ht *hx + 4*c*ht*hx - 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 -pol := lam + lam *(4*s1 *cap + 4*s2 *cap - 3) + lam - - 2 2 2 2 2 2 2 2 2 - *(12*s1 *s2 *cap + 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*lam*cap - - 2 2 2 2 2 2 2 2 2 - *(3*s1 *s2 *cap - 3*s1 *s2 *cap - s1 *cap + s1 - s2 *cap + s2 ) + 8*( - - 2 2 3 2 2 2 2 2 2 2 2 - - 2*s1 *s2 *cap + 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 - - 2 2 2 2 2 2 2 2 2 - (3) 16*cap*(3*s1 *s2 *cap - 3*s1 *s2 *cap - s1 *cap + s1 - s2 *cap + s2 ) - - > 0 - - 2 2 3 2 2 2 2 2 2 2 2 - (4) 8*( - 2*s1 *s2 *cap + 3*s1 *s2 *cap + s1 *cap - 2*s1 *cap + s2 *cap - - 2 - - 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 4 4 - + 8*s1 *s2 *cap - 10*s1 *s2 *cap + 3*s1 *s2 - s1 *cap + s1 - - 2 4 2 2 4 2 4 2 2 - + 8*s1 *s2 *cap - 10*s1 *s2 *cap + 3*s1 *s2 - 2*s1 *s2 *cap - - 2 2 4 4 - + s1 *s2 - s2 *cap + 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 6 4 4 - - 27*s1 *s2 *cap - 32*s1 *s2 *cap + 100*s1 *s2 *cap - - 6 4 3 6 4 2 6 2 4 - - 93*s1 *s2 *cap + 27*s1 *s2 *cap + 10*s1 *s2 *cap - - 6 2 3 6 2 2 6 2 6 3 - - 31*s1 *s2 *cap + 26*s1 *s2 *cap - 6*s1 *s2 *cap - s1 *cap - - 6 2 6 4 6 5 4 6 4 - + 3*s1 *cap - 2*s1 *cap - 32*s1 *s2 *cap + 100*s1 *s2 *cap - - 4 6 3 4 6 2 4 4 4 - - 93*s1 *s2 *cap + 27*s1 *s2 *cap + 20*s1 *s2 *cap - - 4 4 3 4 4 2 4 4 - - 76*s1 *s2 *cap + 73*s1 *s2 *cap - 21*s1 *s2 *cap - - 4 2 3 4 2 2 4 2 - - 3*s1 *s2 *cap + 16*s1 *s2 *cap - 14*s1 *s2 *cap - - 4 2 4 4 2 6 4 - + 3*s1 *s2 - s1 *cap + s1 + 10*s1 *s2 *cap - - 2 6 3 2 6 2 2 6 - - 31*s1 *s2 *cap + 26*s1 *s2 *cap - 6*s1 *s2 *cap - - 2 4 3 2 4 2 2 4 - - 3*s1 *s2 *cap + 16*s1 *s2 *cap - 14*s1 *s2 *cap - - 2 4 2 2 2 2 6 3 6 2 - + 3*s1 *s2 - 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 sin(ax)*c*ht*i -bt := mat((-------------------------------------------,----------------, - hx hx - - sin(ay)*c*ht*i - ----------------), - hx - - sin(ax)*c*ht*i cos(ax)*c*ht - c*ht + hx - (----------------,--------------------------,0), - hx hx - - sin(ay)*c*ht*i 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*ht*hx - - 2 2 2 2 2 2 - - 4*cos(ay)*c *ht + 2*cos(ay)*c*ht*hx + 6*c *ht - 4*c*ht*hx + hx )/hx , - - 2 2 2 2 - sin(ax)*c *ht *i*( - cos(ay) + 1) sin(ay)*c *ht *i*( - cos(ax) + 1) - -----------------------------------,-----------------------------------), - 2 2 - hx hx - - 2 2 - sin(ax)*c *ht *i*(cos(ay) - 1) - (--------------------------------, - 2 - hx - - 2 2 2 2 2 - - 2*cos(ax)*c *ht + 2*cos(ax)*c*ht*hx + 2*c *ht - 2*c*ht*hx + hx - ----------------------------------------------------------------------, - 2 - hx - - 2 2 - sin(ax)*sin(ay)*c *ht - ------------------------), - 2 - hx - - 2 2 2 2 - sin(ay)*c *ht *i*(cos(ax) - 1) sin(ax)*sin(ay)*c *ht - (--------------------------------,------------------------, - 2 2 - hx hx - - 2 2 2 2 2 - - 2*cos(ay)*c *ht + 2*cos(ay)*c*ht*hx + 2*c *ht - 2*c*ht*hx + hx - ----------------------------------------------------------------------)) - 2 - hx - - -bt*b-b*bt; - - - 2 2 - 2*sin(ax)*c *ht *i*( - cos(ay) + 1) -mat((0,-------------------------------------, - 2 - hx - - 2 2 - 2*sin(ay)*c *ht *i*( - cos(ax) + 1) - -------------------------------------), - 2 - hx - - 2 2 - 2*sin(ax)*c *ht *i*(cos(ay) - 1) - (----------------------------------,0,0), - 2 - hx - - 2 2 - 2*sin(ay)*c *ht *i*(cos(ax) - 1) - (----------------------------------,0,0)) - 2 - hx - - -clear a,b,bt,pol; - - - -%*********************************************************************** - -end; -(TIME: fide 15030 15900) +REDUCE 3.6, 15-Jul-95, patched to 6 Mar 96 ... + + +%*********************************************************************** +%***** ***** +%***** Package F I D E - Test Examples Ver. 1.1.2 May 29,1995 ***** +%***** ***** +%*********************************************************************** +%*********************************************************************** +%***** ***** +%***** 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; + + +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) cos(th)*u(2) + 2*sin(th)*u(1) +a2:=------------- + ------------- + df(u(1),r) + ------------------------------- + sin(th)*r r sin(th)*r + + +a3:=curl v; + + df(v(3),th) - df(v(2),fi) cos(th)*v(3) +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) df(v(1),fi,2) +a4:= ( ------------------ + ------------------ + --------------- + df(v(1),r,2) + 2 2 2 2 + sin(th)*r r sin(th) *r + + 2*df(v(1),r) df(v(1),th,2) df(v(1),th)*cos(th) + + -------------- + --------------- + --------------------- + r 2 2 + r sin(th)*r + + 2*(cos(th)*v(2) + sin(th)*v(1)) + - --------------------------------- , + 2 + sin(th)*r + + - 2*df(v(3),fi)*cos(th) df(v(2),fi,2) + -------------------------- + --------------- + df(v(2),r,2) + 2 2 2 2 + sin(th) *r 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 + + 2*df(v(1),th) - v(2) + + --------------- + ------------- , + 2 2 2 + r sin(th) *r + + df(v(3),fi,2) 2*df(v(3),r) df(v(3),th,2) + --------------- + df(v(3),r,2) + -------------- + --------------- + 2 2 r 2 + sin(th) *r r + + df(v(3),th)*cos(th) 2*df(v(2),fi)*cos(th) 2*df(v(1),fi) + + --------------------- + ----------------------- + --------------- + 2 2 2 2 + sin(th)*r sin(th) *r sin(th)*r + + - v(3) + + ------------- ) + 2 2 + sin(th) *r + + +a3:=2*a3+a4; + + - 2*df(v(3),fi) 2*df(v(3),th) - 2*df(v(2),fi) +a3:= ( ------------------ + --------------- + ------------------ + 2 r sin(th)*r + sin(th)*r + + - 2*df(v(2),th) df(v(1),fi,2) 2*df(v(1),r) + + ------------------ + --------------- + df(v(1),r,2) + -------------- + 2 2 2 r + r sin(th) *r + + df(v(1),th,2) df(v(1),th)*cos(th) + + --------------- + --------------------- + 2 2 + r sin(th)*r + + 2*(cos(th)*v(3)*r - cos(th)*v(2) - sin(th)*v(1)) + + -------------------------------------------------- , + 2 + sin(th)*r + + - 2*df(v(3),fi)*cos(th) df(v(2),fi,2) + -------------------------- - 2*df(v(3),r) + --------------- + 2 2 2 2 + sin(th) *r sin(th) *r + + 2*df(v(2),r) df(v(2),th,2) + + df(v(2),r,2) + -------------- + --------------- + r 2 + r + + df(v(2),th)*cos(th) 2*df(v(1),fi) 2*df(v(1),th) + + --------------------- + --------------- + --------------- + 2 sin(th)*r 2 + sin(th)*r r + + 2 + - 2*sin(th) *v(3)*r - v(2) + + ----------------------------- , + 2 2 + sin(th) *r + + df(v(3),fi,2) 2*df(v(3),r) df(v(3),th,2) + --------------- + df(v(3),r,2) + -------------- + --------------- + 2 2 r 2 + sin(th) *r r + + df(v(3),th)*cos(th) 2*df(v(2),fi)*cos(th) + + --------------------- + ----------------------- + 2*df(v(2),r) + 2 2 2 + sin(th)*r sin(th) *r + + 2 + 2*df(v(1),fi) - 2*df(v(1),th) 2*sin(th) *v(2)*r - v(3) + + --------------- + ------------------ + -------------------------- ) + 2 r 2 2 + sin(th)*r sin(th) *r + + +a5:=lapl f; + + df(f,fi,2) 2*df(f,r) df(f,th,2) df(f,th)*cos(th) +a5:=------------- + df(f,r,2) + ----------- + ------------ + ------------------ + 2 2 r 2 2 + sin(th) *r r 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 + + cos(th)*w(2,1) - sin(th)*w(3,3) - sin(th)*w(2,2) + 2*sin(th)*w(1,1) + + --------------------------------------------------------------------- + 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 + + - cos(th)*w(3,3) + cos(th)*w(2,2) + sin(th)*w(2,1) + 2*sin(th)*w(1,2) + ------------------------------------------------------------------------ + 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 + + cos(th)*w(3,2) + cos(th)*w(2,3) + sin(th)*w(3,1) + 2*sin(th)*w(1,3) + + --------------------------------------------------------------------- + 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),fi)*u(3) df(v(1),th)*u(2) + ( ------------------ + df(v(1),r)*u(1) + ------------------ + sin(th)*r r + + u(3)*v(3) + u(2)*v(2) + - ----------------------- , + r + + df(v(2),fi)*u(3) df(v(2),th)*u(2) + ------------------ + df(v(2),r)*u(1) + ------------------ + sin(th)*r r + + - cos(th)*u(3)*v(3) + sin(th)*u(2)*v(1) + + ------------------------------------------ , + sin(th)*r + + df(v(3),fi)*u(3) df(v(3),th)*u(2) + ------------------ + df(v(3),r)*u(1) + ------------------ + sin(th)*r r + + u(3)*(cos(th)*v(2) + sin(th)*v(1)) + + ------------------------------------ ) + 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; + + + +%*********************************************************************** +%***** ***** +%***** 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.2 +***************************** + + Partial Differential Equations + ============================== + +diff(eta,t) - diff(v,x)*eta = 0 + + diff(p,x)*eta +--------------- + diff(v,t) = 0 + ro + + diff(v,x)*eta*p +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 : + ========================================== + +(4*(eta(j,m + 1) - eta(j,m) - eta(j + 1,m) + eta(j + 1,m + 1))*hx - ( + + (eta(j + 1,m + 1) + eta(j,m + 1))*(v(j + 1,m + 1) - v(j,m + 1)) + + + (eta(j + 1,m) + eta(j,m))*(v(j + 1,m) - v(j,m)))*(ht(m + 1) + ht(m)))/(4 + + *(ht(m + 1) + ht(m))*hx) = 0 + + +(4*(v(j,m + 1) - v(j,m) - v(j + 1,m) + v(j + 1,m + 1))*hx*ro + ( + + (eta(j + 1,m + 1) + eta(j,m + 1))*(p(j + 1,m + 1) - p(j,m + 1)) + + + (eta(j + 1,m) + eta(j,m))*(p(j + 1,m) - p(j,m)))*(ht(m + 1) + ht(m)))/(4 + + *(ht(m + 1) + ht(m))*hx*ro) = 0 + + +(4*(eps(j,m + 1) - eps(j,m) - eps(j + 1,m) + eps(j + 1,m + 1))*hx*ro + ( + + (eta(j + 1,m + 1)*p(j + 1,m + 1) + eta(j,m + 1)*p(j,m + 1)) + + *(v(j + 1,m + 1) - v(j,m + 1)) + + + (eta(j + 1,m)*p(j + 1,m) + eta(j,m)*p(j,m))*(v(j + 1,m) - v(j,m))) + + *(ht(m + 1) + ht(m)))/(4*(ht(m + 1) + ht(m))*hx*ro) = 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.2 +***************************** + + 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 : + ========================================== + + 2*j + 1 2*j + 1 2*j + 3 + - ((2*(v(---------,m + 1) + v(---------,m)) - v(---------,m) + 2 2 2 + + 2*j + 3 2*j - 1 2*j - 1 + - v(---------,m + 1) - v(---------,m) - v(---------,m + 1)) + 2 2 2 + + 2*j + 1 + *cos(x(---------))*ht + ( + 2 + + (u(j,m + 1) + u(j,m) - u(j + 1,m) - u(j + 1,m + 1))*ht + + 2 + - (u(j,m + 1) - u(j,m) - u(j + 1,m) + u(j + 1,m + 1))*hx)*hx)/(2*ht*hx ) + + = 0 + + +( - ((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 2 + + + (u(j + 1,m + 1) + u(j + 1,m))*sin(x(j + 1))*ht)/(2*ht*hx) = 0 + + + +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.2 +***************************** + + Partial Differential Equations + ============================== + + diff(ur,x,2) 2 2 + - diff(ui,t) + -------------- = - ur*(ui + ur ) + 2 + + diff(ui,x,2) 2 2 +-------------- + diff(ur,t) = - ui*(ui + ur ) + 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,j + 1) - 2*ur(m,j) + ur(m,j - 1) - 2*ur(m + 1,j) + ur(m + 1,j + 1) + + 2 2 + + ur(m + 1,j - 1))*ht - 4*(ui(m + 1,j) - ui(m,j))*hx )/(4*ht*hx ) = + + 3 3 2 2 + ur(m + 1,j) + ur(m,j) + ui(m,j) *ur(m,j) + ui(m + 1,j) *ur(m + 1,j) + - ----------------------------------------------------------------------- + 2 + + +((ui(m,j + 1) - 2*ui(m,j) + ui(m,j - 1) - 2*ui(m + 1,j) + ui(m + 1,j + 1) + + 2 2 + + ui(m + 1,j - 1))*ht + 4*(ur(m + 1,j) - ur(m,j))*hx )/(4*ht*hx ) = + + 2 2 3 2 + (ui(m,j) + ur(m,j) )*ui(m,j) + ui(m + 1,j) + ui(m + 1,j)*ur(m + 1,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; + + +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.2 +***************************** + + 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 5*diff(u1,x)*p +------------- + ---------------- + ---------------- + ---------------- + 2 2 2 2 + + 5*diff(u2,y)*p + + ---------------- = 0 + 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 *(n(j + 1,i + 1,m + 1)*u1(j + 1,i + 1,m + 1)*hy + + + n(j + 1,i + 1,m + 1)*u2(j + 1,i + 1,m + 1)*hx + + + n(j + 1,i + 1,m)*u1(j + 1,i + 1,m)*hy + + + n(j + 1,i + 1,m)*u2(j + 1,i + 1,m)*hx + + + n(j + 1,i,m + 1)*u1(j + 1,i,m + 1)*hy + + - n(j + 1,i,m + 1)*u2(j + 1,i,m + 1)*hx + + + n(j + 1,i,m)*u1(j + 1,i,m)*hy - n(j + 1,i,m)*u2(j + 1,i,m)*hx + + - n(j,i + 1,m + 1)*u1(j,i + 1,m + 1)*hy + + + n(j,i + 1,m + 1)*u2(j,i + 1,m + 1)*hx + + - n(j,i + 1,m)*u1(j,i + 1,m)*hy + n(j,i + 1,m)*u2(j,i + 1,m)*hx + + - n(j,i,m + 1)*u1(j,i,m + 1)*hy - n(j,i,m + 1)*u2(j,i,m + 1)*hx + + -1 + - n(j,i,m)*u1(j,i,m)*hy - n(j,i,m)*u2(j,i,m)*hx))/4 + (ht *( + + n(j,i,m + 1) - n(j,i,m) - n(j,i + 1,m) + n(j,i + 1,m + 1) - n(j + 1,i,m) + + + n(j + 1,i,m + 1) - n(j + 1,i + 1,m) + n(j + 1,i + 1,m + 1)))/4 = 0 + + + -1 +(hx *((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)) + + + (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 + 1,m + 1)*u1(j + 1,i + 1,m + 1) + + + n(j,i + 1,m + 1)*u1(j,i + 1,m + 1)) + + -1 -1 -1 + *(u1(j + 1,i + 1,m + 1) - u1(j,i + 1,m + 1)))*m)/8 + (hy *hx *ht *((( + + (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)) + + + (n(j + 1,i,m + 1) + n(j + 1,i,m)) + + *(u1(j + 1,i,m + 1) - u1(j + 1,i,m)) + + + (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)))*hx*m + 2*( + + p(j + 1,i,m + 1) + p(j + 1,i,m) + p(j + 1,i + 1,m) + + + p(j + 1,i + 1,m + 1) + + - (p(j,i,m + 1) + p(j,i,m) + p(j,i + 1,m) + p(j,i + 1,m + 1)))*ht) + + *hy + ((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)) + + *(u1(j,i + 1,m) - u1(j,i,m)) + + + (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 + 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)))*ht*hx*m))/8 = 0 + + + -1 -1 +(hy *hx *(((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)) + + *(u2(j + 1,i,m) - u2(j,i,m)) + + + (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 + 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)))*hy + ( + + (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)) + + *(u2(j,i + 1,m) - u2(j,i,m)) + + + (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 + 1,i + 1,m + 1)*u2(j + 1,i + 1,m + 1) + + + n(j + 1,i,m + 1)*u2(j + 1,i,m + 1)) + + -1 + *(u2(j + 1,i + 1,m + 1) - u2(j + 1,i,m + 1)))*hx)*m)/8 + ( - hy + + -1 + *ht *(2*(p(j,i,m + 1) + p(j,i,m) - p(j,i + 1,m) - p(j,i + 1,m + 1) + + + p(j + 1,i,m) + p(j + 1,i,m + 1) - p(j + 1,i + 1,m) + + - p(j + 1,i + 1,m + 1))*ht - ((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)) + + + (n(j + 1,i,m + 1) + n(j + 1,i,m)) + + *(u2(j + 1,i,m + 1) - u2(j + 1,i,m)) + + + (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)))*hy*m))/8 = 0 + + + -1 -1 +(hy *hx *(3*((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)) + + + (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 + 1,m + 1) - p(j,i + 1,m + 1)) + + *(u1(j + 1,i + 1,m + 1) + u1(j,i + 1,m + 1)))*hy + 2*( + + 4*p(j + 1,i + 1,m + 1)*u2(j + 1,i + 1,m + 1) + + - p(j + 1,i + 1,m + 1)*u2(j + 1,i,m + 1) + + + 4*p(j + 1,i + 1,m)*u2(j + 1,i + 1,m) + + - p(j + 1,i + 1,m)*u2(j + 1,i,m) + + + p(j + 1,i,m + 1)*u2(j + 1,i + 1,m + 1) + + - 4*p(j + 1,i,m + 1)*u2(j + 1,i,m + 1) + + + p(j + 1,i,m)*u2(j + 1,i + 1,m) - 4*p(j + 1,i,m)*u2(j + 1,i,m) + + + 4*p(j,i + 1,m + 1)*u2(j,i + 1,m + 1) + + - p(j,i + 1,m + 1)*u2(j,i,m + 1) + 4*p(j,i + 1,m)*u2(j,i + 1,m) + + - p(j,i + 1,m)*u2(j,i,m) + p(j,i,m + 1)*u2(j,i + 1,m + 1) + + - 4*p(j,i,m + 1)*u2(j,i,m + 1) + p(j,i,m)*u2(j,i + 1,m) + + - 4*p(j,i,m)*u2(j,i,m))*hx + 5*( + + (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)) + + + (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 + 1,m + 1) + p(j,i + 1,m + 1)) + + -1 + *(u1(j + 1,i + 1,m + 1) - u1(j,i + 1,m + 1)))*hy))/16 + (3*ht *( + + p(j,i,m + 1) - p(j,i,m) - p(j,i + 1,m) + p(j,i + 1,m + 1) - p(j + 1,i,m) + + + p(j + 1,i,m + 1) - p(j + 1,i + 1,m) + p(j + 1,i + 1,m + 1)))/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.2 +***************************** + + 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 + + 3*diff(tt,t)*k*n 3*diff(tt,x)*k*n*u +diff(q,x) + ------------------ + -------------------- + 2 2 + + diff(u,x)*(2*k*n*tt + p) + + -------------------------- = 0 + 2 + + 2*diff(q,x) +diff(p,t) + diff(p,x)*u + ------------- + diff(u,x)*(k*n*tt + p) = 0 + 5 + + 2 + 5*diff(tt,x)*k *n*tt +diff(p,x)*(k*n*tt - p) + diff(q,t) + diff(q,x)*u + ---------------------- + 2*m + + + diff(u,x)*q = 0 + + +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) - (u(j,m + 1) + u(j,m)))))/4 + + -1 + ht *(n(j,m + 1) - n(j,m) - n(j + 1,m) + n(j + 1,m + 1)) + + ---------------------------------------------------------- = 0 + 2 + + + -1 +( - hx *(n(j,m + 1)*tt(j,m + 1) + n(j,m)*tt(j,m) - n(j + 1,m)*tt(j + 1,m) + + -1 -1 + - n(j + 1,m + 1)*tt(j + 1,m + 1))*k)/2 + (hx *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)))*hx*m + + + 2*(p(j + 1,m + 1) + p(j + 1,m) - (p(j,m + 1) + p(j,m)))*ht + ( + + (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)))*ht*m) + + )/4 = 0 + + + -1 +(hx *((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)/4 + + -1 -1 ++ (hx *ht *(((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) - (q(j,m + 1) + q(j,m))))*ht + + + 3*((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)))*hx*k + 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)) + + )*ht*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 + + + 2*(5*p(j + 1,m + 1)*u(j + 1,m + 1) + 5*p(j + 1,m)*u(j + 1,m) + + - 5*p(j,m + 1)*u(j,m + 1) - 5*p(j,m)*u(j,m) + 2*q(j + 1,m + 1) + + + 2*q(j + 1,m) - 2*q(j,m + 1) - 2*q(j,m))))/20 + + -1 + ht *(p(j,m + 1) - p(j,m) - p(j + 1,m) + p(j + 1,m + 1)) + + ---------------------------------------------------------- = 0 + 2 + + + -1 +( - hx *(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)) - 2*( + + q(j + 1,m + 1)*u(j + 1,m + 1) + q(j + 1,m)*u(j + 1,m) + + - q(j,m + 1)*u(j,m + 1) - q(j,m)*u(j,m)))*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))*(tt(j + 1,m) - tt(j,m))) + + 2 + *k - 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))/(8*m) + + -1 + ht *(q(j,m + 1) - q(j,m) - q(j + 1,m) + q(j + 1,m + 1)) + + ---------------------------------------------------------- = 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; + +maxorder x=4; + + +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: + + 2 +hx + + 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] + + + [ hx 2*sin(ax)*ht*( - i*cos(ax) + sin(ax))] +h0 := [ ] + [ - 2*i*sin(ax)*ht hx ] + + + [ 2*sin(ax)*ht*( - i*cos(ax) + sin(ax)) ] + [ 1 ---------------------------------------] + [ hx ] +bb := [ ] + [ - 2*i*sin(ax)*ht ] + [------------------- 1 ] + [ hx ] + + +bb:=denotemat bb; + + + [ 1 ai12*i + ar12] +bb := [ ] + [ai21*i 1 ] + + +factor lam; + + +pol:=charpol bb; + + + 2 +pol := lam - 2*lam + ai12*ai21 - i*ai21*ar12 + 1 + +prdenot; + + + 2 + 2*sin(ax) *ht +ar12 := --------------- + hx + + - 2*cos(ax)*sin(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 := [cos(ax)*a*ht - i*sin(ax)*a*ht - a*ht + hx] + + + [ cos(ax)*a*ht - i*sin(ax)*a*ht - a*ht + hx ] +bb := [-------------------------------------------] + [ hx ] + + +bb:=denotemat bb; + + +bb := [ai11*i + 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 + + + [ i*sin(ax)*ht*r0 i*sin(ax)*ht*u0 + hx] +h1 := [ ] + [2*r0*(i*sin(ax)*ht*u0 + hx) 0 ] + + +h0 := mat((0,hx), + + 2 2 + (cos(ax)*c0 *ht - i*sin(ax)*c0 *ht + 2*hx*r0, + + 2 + c0 *ht*( - cos(ax) - i*sin(ax)))) + + + 2 2 + - i*cos(ax)*c0 *ht - sin(ax)*c0 *ht - 2*i*hx*r0 +bb := mat((--------------------------------------------------, + 2*r0*(sin(ax)*ht*u0 - i*hx) + + 2 + c0 *ht*(i*cos(ax) - sin(ax)) + ------------------------------), + 2*r0*(sin(ax)*ht*u0 - i*hx) + + 2 2 + sin(ax)*ht*(i*cos(ax)*c0 *ht + sin(ax)*c0 *ht + 2*i*hx*r0) + (------------------------------------------------------------,( + 2 2 2 2 + 2*(sin(ax) *ht *u0 - 2*i*sin(ax)*ht*hx*u0 - hx ) + + 2 2 2 2 2 + - i*cos(ax)*sin(ax)*c0 *ht + sin(ax) *c0 *ht + + 2 + - 2*i*sin(ax)*ht*hx*u0 - 2*hx )/(2 + + 2 2 2 2 + *(sin(ax) *ht *u0 - 2*i*sin(ax)*ht*hx*u0 - hx )))) + + +bb:=denotemat bb; + + + [ai11*i + ar11 ai12*i + ar12] +bb := [ ] + [ai21*i + ar21 ai22*i + ar22] + + +pol:=charpol bb; + + + 2 +pol := lam + lam*( - i*ai11 - i*ai22 - ar11 - ar22) - ai11*ai22 + i*ai11*ar22 + + + ai12*ai21 - i*ai12*ar21 - i*ai21*ar12 + i*ai22*ar11 + ar11*ar22 + + - ar12*ar21 + +prdenot; + + + 2 + - c0 +ar11 := --------- + 2*r0*u0 + + 2 2 + - cos(ax)*c0 *ht*u0 - c0 *hx - 2*hx*r0*u0 +ai11 := -------------------------------------------- + 2*r0*u0*(sin(ax)*ht*u0 - hx*i) + + 2 + - c0 +ar12 := --------- + 2*r0*u0 + + 2 + c0 *(cos(ax)*ht*u0 - hx) +ai12 := -------------------------------- + 2*r0*u0*(sin(ax)*ht*u0 - hx*i) + + 2 2 2 + sin(ax) *c0 *ht +ar21 := ---------------------------- + 2 2 2 2 + 2*(sin(ax) *ht *u0 - hx ) + + 2 2 3 2 2 2 +ai21 := (sin(ax)*ht*(cos(ax)*sin(ax) *c0 *ht *u0 - cos(ax)*c0 *ht*hx + + 2 2 2 2 2 2 3 + + 2*sin(ax) *c0 *ht *hx*u0 + 2*sin(ax) *ht *hx*r0*u0 - 2*hx *r0))/ + + 4 4 4 3 3 3 2 2 2 2 + (2*(sin(ax) *ht *u0 - 2*sin(ax) *ht *hx*i*u0 - 2*sin(ax) *ht *hx *u0 + + 3 4 + + 2*sin(ax)*ht*hx *i*u0 + hx )) + + 2 2 2 2 + sin(ax) *c0 *ht - 2*hx +ar22 := ---------------------------- + 2 2 2 2 + 2*(sin(ax) *ht *u0 - hx ) + + 2 2 3 2 2 2 +ai22 := (sin(ax)*ht*( - cos(ax)*sin(ax) *c0 *ht *u0 + cos(ax)*c0 *ht*hx + + 2 2 2 2 2 3 3 + + 2*sin(ax) *c0 *ht *hx*u0 - 2*sin(ax) *ht *hx*u0 - 2*hx *u0))/(2* + + 4 4 4 3 3 3 2 2 2 2 + (sin(ax) *ht *u0 - 2*sin(ax) *ht *hx*i*u0 - 2*sin(ax) *ht *hx *u0 + + 3 4 + + 2*sin(ax)*ht*hx *i*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 + + + [ hx - i*sin(ax)*c*ht ] +h1 := [ ] + [sin(ax)*c*ht*( - i*cos(ax) - sin(ax)) hx*(cos(ax) - i*sin(ax))] + + + [ hx i*sin(ax)*c*ht ] +h0 := [ ] + [sin(ax)*c*ht*(i*cos(ax) + sin(ax)) hx*(cos(ax) - i*sin(ax))] + + + [ 2 2 2 2 ] + [ - sin(ax) *c *ht + hx 2*i*sin(ax)*c*ht*hx ] + [-------------------------- ----------------------- ] + [ 2 2 2 2 2 2 2 2 ] + [ sin(ax) *c *ht + hx sin(ax) *c *ht + hx ] +bb := [ ] + [ 2 2 2 2 ] + [ 2*i*sin(ax)*c*ht*hx - 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 ai12*i] +bb := [ ] + [ai21*i 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*ht*hx +ai12 := ----------------------- + 2 2 2 2 + sin(ax) *c *ht + hx + + 2*sin(ax)*c*ht*hx +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 + + + cos(ax)*c*ht*hy + cos(ay)*c*ht*hx - c*ht*hx - c*ht*hy + hx*hy +bb := mat((---------------------------------------------------------------, + hx*hy + + - i*sin(ax)*c*ht - i*sin(ay)*c*ht + -------------------,-------------------), + hx hy + + - i*sin(ax)*c*ht cos(ax)*c*ht - c*ht + hx + (-------------------,--------------------------,0), + hx hx + + - i*sin(ay)*c*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*ht*hx + + + 2*c*ht*hx + 2*c*ht*hy - 3*hx*hy) + lam*( + + 2 2 2 2 + 3*cos(ax)*cos(ay)*c *ht *hx*hy - 3*cos(ax)*c *ht *hx*hy + + 2 2 2 2 2 2 2 2 + - 2*cos(ax)*c *ht *hy + 4*cos(ax)*c*ht*hx*hy + cos(ay) *c *ht *hx + + 2 2 2 2 2 + - 2*cos(ay)*c *ht *hx - 3*cos(ay)*c *ht *hx*hy + + 2 2 2 2 2 2 2 2 + + 4*cos(ay)*c*ht*hx *hy + sin(ay) *c *ht *hx + c *ht *hx + + 2 2 2 2 2 2 2 + + 3*c *ht *hx*hy + 2*c *ht *hy - 4*c*ht*hx *hy - 4*c*ht*hx*hy + + 2 2 2 3 3 + + 3*hx *hy ) - cos(ax)*cos(ay) *c *ht *hx + + 3 3 3 3 + + 2*cos(ax)*cos(ay)*c *ht *hx + 2*cos(ax)*cos(ay)*c *ht *hy + + 2 2 2 3 3 + - 3*cos(ax)*cos(ay)*c *ht *hx*hy - cos(ax)*sin(ay) *c *ht *hx + + 3 3 3 3 2 2 + - cos(ax)*c *ht *hx - 2*cos(ax)*c *ht *hy + 3*cos(ax)*c *ht *hx*hy + + 2 2 2 2 2 3 3 + + 2*cos(ax)*c *ht *hy - 2*cos(ax)*c*ht*hx*hy + cos(ay) *c *ht *hx + + 2 2 2 2 3 3 3 3 + - cos(ay) *c *ht *hx - 2*cos(ay)*c *ht *hx - 2*cos(ay)*c *ht *hy + + 2 2 2 2 2 2 + + 2*cos(ay)*c *ht *hx + 3*cos(ay)*c *ht *hx*hy - 2*cos(ay)*c*ht*hx *hy + + 2 3 3 2 2 2 2 3 3 3 3 + + sin(ay) *c *ht *hx - sin(ay) *c *ht *hx + c *ht *hx + 2*c *ht *hy + + 2 2 2 2 2 2 2 2 2 + - c *ht *hx - 3*c *ht *hx*hy - 2*c *ht *hy + 2*c*ht*hx *hy + + 2 2 2 2 2 + + 2*c*ht*hx*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 2 2 +pol := lam + lam *(4*s1 *cap1 + 4*s2 *cap2 - 3) + lam*(12*s1 *s2 *cap1*cap2 + + 2 2 2 2 2 2 + + 4*s1 *cap1 - 8*s1 *cap1 + 4*s2 *cap2 - 8*s2 *cap2 + 3) + + 2 2 2 2 2 2 2 2 + + 8*s1 *s2 *cap1 *cap2 + 8*s1 *s2 *cap1*cap2 - 12*s1 *s2 *cap1*cap2 + + 2 2 2 2 2 2 + - 4*s1 *cap1 + 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 2 2 +pol := lam + lam *(4*s1 *cap1 + 4*s2 *cap2 - 3) + lam*(12*s1 *s2 *cap1*cap2 + + 2 2 2 2 2 2 + + 4*s1 *cap1 - 8*s1 *cap1 + 4*s2 *cap2 - 8*s2 *cap2 + 3) + + 2 2 2 2 2 2 2 2 + + 8*s1 *s2 *cap1 *cap2 + 8*s1 *s2 *cap1*cap2 - 12*s1 *s2 *cap1*cap2 + + 2 2 2 2 2 2 + - 4*s1 *cap1 + 4*s1 *cap1 - 4*s2 *cap2 + 4*s2 *cap2 - 1 + +pol1:=hurw pol; + + + 3 2 2 2 2 2 2 +pol1 := 8*lam *s1 *s2 *cap1*cap2*(cap1 + cap2) + 8*lam *( - 3*s1 *s2 *cap1 *cap2 + + 2 2 2 2 2 2 2 2 2 + - 3*s1 *s2 *cap1*cap2 + 3*s1 *s2 *cap1*cap2 + s1 *cap1 + s2 *cap2 + + 2 2 2 2 2 2 + ) + 8*lam*(3*s1 *s2 *cap1 *cap2 + 3*s1 *s2 *cap1*cap2 + + 2 2 2 2 2 2 2 + - 6*s1 *s2 *cap1*cap2 - 2*s1 *cap1 + 2*s1 *cap1 - 2*s2 *cap2 + + 2 2 2 2 2 2 2 + + 2*s2 *cap2) + 8*( - s1 *s2 *cap1 *cap2 - s1 *s2 *cap1*cap2 + + 2 2 2 2 2 2 2 + + 3*s1 *s2 *cap1*cap2 + s1 *cap1 - 2*s1 *cap1 + s2 *cap2 + + 2 + - 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 2 2 +cpr00 := 8*s1 *s2 *cap1 *cap2 + 8*s1 *s2 *cap1*cap2 - 12*s1 *s2 *cap1*cap2 + + 2 2 2 2 2 2 + - 4*s1 *cap1 + 4*s1 *cap1 - 4*s2 *cap2 + 4*s2 *cap2 - 1 + +cpr01 := + + 2 2 2 2 2 2 2 2 +12*s1 *s2 *cap1*cap2 + 4*s1 *cap1 - 8*s1 *cap1 + 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 + + +bb := mat((( - i*cos(ax)*sin(ax)*ht*hx*w - i*cos(ax)*sin(ay)*ht*hx*v + + - i*cos(ay)*sin(ax)*ht*hx*w - i*cos(ay)*sin(ay)*ht*hx*v + + 2 2 2 2 2 2 2 + - 2*sin(ax) *ht *w - 4*sin(ax)*sin(ay)*ht *v*w - 2*sin(ay) *ht *v + + 2 2 + + hx )/hx ,(sin(ax)*ht*p*( - i*cos(ax)*hx - i*cos(ay)*hx + + 2 + - 4*sin(ax)*ht*w - 4*sin(ay)*ht*v))/hx ,(ht*p*( + + - i*cos(ax)*sin(ay)*hx - 4*i*cos(ay)*sin(ay)*ht*v + + 2 + - i*cos(ay)*sin(ay)*hx - 4*sin(ax)*sin(ay)*ht*w - 2*ht*v))/hx + + 2 2 2 + - 2*ht *(sin(ax) + sin(ay) ) + ,--------------------------------), + 2 + hx + + (0,( - i*cos(ax)*sin(ax)*ht*hx*w - i*cos(ax)*sin(ay)*ht*hx*v + + - i*cos(ay)*sin(ax)*ht*hx*w - i*cos(ay)*sin(ay)*ht*hx*v + + 2 2 2 2 2 2 + - 2*sin(ax) *c *ht - 2*sin(ax) *ht *w + + 2 2 2 2 2 2 + - 4*sin(ax)*sin(ay)*ht *v*w - 2*sin(ay) *ht *v + hx )/hx , + + 2 2 + - 2*sin(ax)*sin(ay)*c *ht + -----------------------------,(sin(ax)*ht*( - i*cos(ax)*hx + 2 + hx + + 2 + - i*cos(ay)*hx - 4*sin(ax)*ht*w - 4*sin(ay)*ht*v))/(hx *p)), + + 2 2 + - 2*sin(ax)*sin(ay)*c *ht + (0,-----------------------------,( - i*cos(ax)*sin(ax)*ht*hx*w + 2 + hx + + - i*cos(ax)*sin(ay)*ht*hx*v - i*cos(ay)*sin(ax)*ht*hx*w + + 2 2 2 + - i*cos(ay)*sin(ay)*ht*hx*v - 2*sin(ax) *ht *w + + 2 2 2 2 + - 4*sin(ax)*sin(ay)*ht *v*w - 2*sin(ay) *c *ht + + 2 2 2 2 2 + - 2*sin(ay) *ht *v + hx )/hx ,(ht*( - 2*cos(ax)*cos(ay)*ht*w + + - 2*i*cos(ax)*sin(ay)*ht*w - i*cos(ax)*sin(ay)*hx + + - 2*i*cos(ay)*sin(ax)*ht*w - i*cos(ay)*sin(ay)*hx + + 2 2 + - 2*sin(ax)*sin(ay)*ht*w - 4*sin(ay) *ht*v))/(hx *p)), + + 2 + (0,(sin(ax)*c *ht*p*( - i*cos(ax)*hx - i*cos(ay)*hx - 4*sin(ax)*ht*w + + 2 2 + - 4*sin(ay)*ht*v))/hx ,(sin(ay)*c *ht*p*( - i*cos(ax)*hx + + 2 + - i*cos(ay)*hx - 4*sin(ax)*ht*w - 4*sin(ay)*ht*v))/hx ,( + + - i*cos(ax)*sin(ax)*ht*hx*w - i*cos(ax)*sin(ay)*ht*hx*v + + - i*cos(ay)*sin(ax)*ht*hx*w - i*cos(ay)*sin(ay)*ht*hx*v + + 2 2 2 2 2 2 + - 2*sin(ax) *c *ht - 2*sin(ax) *ht *w + + 2 2 2 2 + - 4*sin(ax)*sin(ay)*ht *v*w - 2*sin(ay) *c *ht + + 2 2 2 2 2 + - 2*sin(ay) *ht *v + 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; + + + [ai11*i + ar11 ai12*i + ar12 ai13*i + ar13 ar14 ] + [ ] + [ 0 ai22*i + ar22 ar23 ai24*i + ar24] +bb := [ ] + [ 0 ar32 ai33*i + ar33 ai34*i + ar34] + [ ] + [ 0 ai42*i + ar42 ai43*i + ar43 ai44*i + 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 - ar33 - ar44) + + + 2 + lam *( - ai11*ai22 - ai11*ai33 - ai11*ai44 + i*ai11*ar22 + i*ai11*ar33 + + + i*ai11*ar44 - ai22*ai33 - ai22*ai44 + i*ai22*ar11 + i*ai22*ar33 + + + i*ai22*ar44 + ai24*ai42 - i*ai24*ar42 - ai33*ai44 + i*ai33*ar11 + + + i*ai33*ar22 + i*ai33*ar44 + ai34*ai43 - i*ai34*ar43 + + - i*ai42*ar24 - i*ai43*ar34 + i*ai44*ar11 + i*ai44*ar22 + + + i*ai44*ar33 + ar11*ar22 + ar11*ar33 + ar11*ar44 + ar22*ar33 + + + ar22*ar44 - ar23*ar32 - ar24*ar42 + ar33*ar44 - ar34*ar43) + lam + + *(i*ai11*ai22*ai33 + i*ai11*ai22*ai44 + ai11*ai22*ar33 + ai11*ai22*ar44 + + - i*ai11*ai24*ai42 - ai11*ai24*ar42 + i*ai11*ai33*ai44 + + + ai11*ai33*ar22 + ai11*ai33*ar44 - i*ai11*ai34*ai43 - ai11*ai34*ar43 + + - ai11*ai42*ar24 - ai11*ai43*ar34 + ai11*ai44*ar22 + ai11*ai44*ar33 + + - i*ai11*ar22*ar33 - i*ai11*ar22*ar44 + i*ai11*ar23*ar32 + + + i*ai11*ar24*ar42 - i*ai11*ar33*ar44 + i*ai11*ar34*ar43 + + + i*ai22*ai33*ai44 + ai22*ai33*ar11 + ai22*ai33*ar44 + + - i*ai22*ai34*ai43 - ai22*ai34*ar43 - ai22*ai43*ar34 + ai22*ai44*ar11 + + + ai22*ai44*ar33 - i*ai22*ar11*ar33 - i*ai22*ar11*ar44 + + - i*ai22*ar33*ar44 + i*ai22*ar34*ar43 - i*ai24*ai33*ai42 + + - ai24*ai33*ar42 - ai24*ai42*ar11 - ai24*ai42*ar33 + ai24*ai43*ar32 + + + i*ai24*ar11*ar42 - i*ai24*ar32*ar43 + i*ai24*ar33*ar42 + + - ai33*ai42*ar24 + ai33*ai44*ar11 + ai33*ai44*ar22 - i*ai33*ar11*ar22 + + - i*ai33*ar11*ar44 - i*ai33*ar22*ar44 + i*ai33*ar24*ar42 + + + ai34*ai42*ar23 - ai34*ai43*ar11 - ai34*ai43*ar22 + i*ai34*ar11*ar43 + + + i*ai34*ar22*ar43 - i*ai34*ar23*ar42 + i*ai42*ar11*ar24 + + - i*ai42*ar23*ar34 + i*ai42*ar24*ar33 + i*ai43*ar11*ar34 + + + i*ai43*ar22*ar34 - i*ai43*ar24*ar32 - i*ai44*ar11*ar22 + + - i*ai44*ar11*ar33 - i*ai44*ar22*ar33 + i*ai44*ar23*ar32 + + - ar11*ar22*ar33 - ar11*ar22*ar44 + ar11*ar23*ar32 + ar11*ar24*ar42 + + - ar11*ar33*ar44 + ar11*ar34*ar43 - ar22*ar33*ar44 + ar22*ar34*ar43 + + + ar23*ar32*ar44 - ar23*ar34*ar42 - ar24*ar32*ar43 + ar24*ar33*ar42) + + + ai11*ai22*ai33*ai44 - i*ai11*ai22*ai33*ar44 - ai11*ai22*ai34*ai43 + + + i*ai11*ai22*ai34*ar43 + i*ai11*ai22*ai43*ar34 - i*ai11*ai22*ai44*ar33 + + - ai11*ai22*ar33*ar44 + ai11*ai22*ar34*ar43 - ai11*ai24*ai33*ai42 + + + i*ai11*ai24*ai33*ar42 + i*ai11*ai24*ai42*ar33 - i*ai11*ai24*ai43*ar32 + + - ai11*ai24*ar32*ar43 + ai11*ai24*ar33*ar42 + i*ai11*ai33*ai42*ar24 + + - i*ai11*ai33*ai44*ar22 - ai11*ai33*ar22*ar44 + ai11*ai33*ar24*ar42 + + - i*ai11*ai34*ai42*ar23 + i*ai11*ai34*ai43*ar22 + ai11*ai34*ar22*ar43 + + - ai11*ai34*ar23*ar42 - ai11*ai42*ar23*ar34 + ai11*ai42*ar24*ar33 + + + ai11*ai43*ar22*ar34 - ai11*ai43*ar24*ar32 - ai11*ai44*ar22*ar33 + + + ai11*ai44*ar23*ar32 + i*ai11*ar22*ar33*ar44 - i*ai11*ar22*ar34*ar43 + + - i*ai11*ar23*ar32*ar44 + i*ai11*ar23*ar34*ar42 + i*ai11*ar24*ar32*ar43 + + - i*ai11*ar24*ar33*ar42 - i*ai22*ai33*ai44*ar11 - ai22*ai33*ar11*ar44 + + + i*ai22*ai34*ai43*ar11 + ai22*ai34*ar11*ar43 + ai22*ai43*ar11*ar34 + + - ai22*ai44*ar11*ar33 + i*ai22*ar11*ar33*ar44 - i*ai22*ar11*ar34*ar43 + + + i*ai24*ai33*ai42*ar11 + ai24*ai33*ar11*ar42 + ai24*ai42*ar11*ar33 + + - ai24*ai43*ar11*ar32 + i*ai24*ar11*ar32*ar43 - i*ai24*ar11*ar33*ar42 + + + ai33*ai42*ar11*ar24 - ai33*ai44*ar11*ar22 + i*ai33*ar11*ar22*ar44 + + - i*ai33*ar11*ar24*ar42 - ai34*ai42*ar11*ar23 + ai34*ai43*ar11*ar22 + + - i*ai34*ar11*ar22*ar43 + i*ai34*ar11*ar23*ar42 + i*ai42*ar11*ar23*ar34 + + - i*ai42*ar11*ar24*ar33 - i*ai43*ar11*ar22*ar34 + i*ai43*ar11*ar24*ar32 + + + i*ai44*ar11*ar22*ar33 - i*ai44*ar11*ar23*ar32 + ar11*ar22*ar33*ar44 + + - ar11*ar22*ar34*ar43 - ar11*ar23*ar32*ar44 + ar11*ar23*ar34*ar42 + + + ar11*ar24*ar32*ar43 - ar11*ar24*ar33*ar42 + +denotid cp; + + +pol:=denotepol pol; + + + 4 3 2 +pol := lam + lam *(cpi03*i + cpr03) + lam *(cpi02*i + cpr02) + + + lam*(cpi01*i + cpr01) + cpi00*i + 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 *(cpi03 + 2*cpr02 + cpr03 ) + + 5 + + 2*lam *(cpi02*cpi03 + cpr01 + cpr02*cpr03) + + 4 2 2 + + lam *(2*cpi01*cpi03 + cpi02 + 2*cpr00 + 2*cpr01*cpr03 + cpr02 ) + + 3 + + 2*lam *(cpi00*cpi03 + cpi01*cpi02 + cpr00*cpr03 + cpr01*cpr02) + + 2 2 2 + + lam *(2*cpi00*cpi02 + cpi01 + 2*cpr00*cpr02 + cpr01 ) + + 2 2 + + 2*lam*(cpi00*cpi01 + cpr00*cpr01) + cpi00 + cpr00 + +denotid rp; + + +pol:=denotepol pol; + + + 8 7 6 5 4 3 +pol := lam + lam *rpr07 + lam *rpr06 + lam *rpr05 + lam *rpr04 + lam *rpr03 + + 2 + + 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*k3 *p*(s1*k1 + s2*k2) +ar42 := ------------------------------- + r1 + + 2 + - s1*k3 *p*(c1 + c2) +ai42 := ----------------------- + r1 + + 2 + - 4*s2*k3 *p*(s1*k1 + s2*k2) +ar43 := ------------------------------- + r1 + + 2 + - s2*k3 *p*(c1 + c2) +ai43 := ----------------------- + r1 + + 2 2 2 2 2 2 2 2 +ar44 := - 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 := ai11*ai22*ai33*ai44 - ai11*ai22*ai34*ai43 - ai11*ai22*ar33*ar44 + + + ai11*ai22*ar34*ar43 - ai11*ai24*ai33*ai42 - ai11*ai24*ar32*ar43 + + + ai11*ai24*ar33*ar42 - ai11*ai33*ar22*ar44 + ai11*ai33*ar24*ar42 + + + ai11*ai34*ar22*ar43 - ai11*ai34*ar23*ar42 - ai11*ai42*ar23*ar34 + + + ai11*ai42*ar24*ar33 + ai11*ai43*ar22*ar34 - ai11*ai43*ar24*ar32 + + - ai11*ai44*ar22*ar33 + ai11*ai44*ar23*ar32 - ai22*ai33*ar11*ar44 + + + ai22*ai34*ar11*ar43 + ai22*ai43*ar11*ar34 - ai22*ai44*ar11*ar33 + + + ai24*ai33*ar11*ar42 + ai24*ai42*ar11*ar33 - ai24*ai43*ar11*ar32 + + + ai33*ai42*ar11*ar24 - ai33*ai44*ar11*ar22 - ai34*ai42*ar11*ar23 + + + ai34*ai43*ar11*ar22 + ar11*ar22*ar33*ar44 - ar11*ar22*ar34*ar43 + + - ar11*ar23*ar32*ar44 + ar11*ar23*ar34*ar42 + ar11*ar24*ar32*ar43 + + - ar11*ar24*ar33*ar42 + +cpi00 := - ai11*ai22*ai33*ar44 + ai11*ai22*ai34*ar43 + ai11*ai22*ai43*ar34 + + - ai11*ai22*ai44*ar33 + ai11*ai24*ai33*ar42 + ai11*ai24*ai42*ar33 + + - ai11*ai24*ai43*ar32 + ai11*ai33*ai42*ar24 - ai11*ai33*ai44*ar22 + + - ai11*ai34*ai42*ar23 + ai11*ai34*ai43*ar22 + ai11*ar22*ar33*ar44 + + - ai11*ar22*ar34*ar43 - ai11*ar23*ar32*ar44 + ai11*ar23*ar34*ar42 + + + ai11*ar24*ar32*ar43 - ai11*ar24*ar33*ar42 - ai22*ai33*ai44*ar11 + + + ai22*ai34*ai43*ar11 + ai22*ar11*ar33*ar44 - ai22*ar11*ar34*ar43 + + + ai24*ai33*ai42*ar11 + ai24*ar11*ar32*ar43 - ai24*ar11*ar33*ar42 + + + ai33*ar11*ar22*ar44 - ai33*ar11*ar24*ar42 - ai34*ar11*ar22*ar43 + + + ai34*ar11*ar23*ar42 + ai42*ar11*ar23*ar34 - ai42*ar11*ar24*ar33 + + - ai43*ar11*ar22*ar34 + ai43*ar11*ar24*ar32 + ai44*ar11*ar22*ar33 + + - ai44*ar11*ar23*ar32 + +cpr01 := ai11*ai22*ar33 + ai11*ai22*ar44 - ai11*ai24*ar42 + ai11*ai33*ar22 + + + ai11*ai33*ar44 - ai11*ai34*ar43 - ai11*ai42*ar24 - ai11*ai43*ar34 + + + ai11*ai44*ar22 + ai11*ai44*ar33 + ai22*ai33*ar11 + ai22*ai33*ar44 + + - ai22*ai34*ar43 - ai22*ai43*ar34 + ai22*ai44*ar11 + ai22*ai44*ar33 + + - ai24*ai33*ar42 - ai24*ai42*ar11 - ai24*ai42*ar33 + ai24*ai43*ar32 + + - ai33*ai42*ar24 + ai33*ai44*ar11 + ai33*ai44*ar22 + ai34*ai42*ar23 + + - ai34*ai43*ar11 - ai34*ai43*ar22 - ar11*ar22*ar33 - ar11*ar22*ar44 + + + ar11*ar23*ar32 + ar11*ar24*ar42 - ar11*ar33*ar44 + ar11*ar34*ar43 + + - ar22*ar33*ar44 + ar22*ar34*ar43 + ar23*ar32*ar44 - ar23*ar34*ar42 + + - ar24*ar32*ar43 + ar24*ar33*ar42 + +cpi01 := ai11*ai22*ai33 + ai11*ai22*ai44 - ai11*ai24*ai42 + ai11*ai33*ai44 + + - ai11*ai34*ai43 - ai11*ar22*ar33 - ai11*ar22*ar44 + ai11*ar23*ar32 + + + ai11*ar24*ar42 - ai11*ar33*ar44 + ai11*ar34*ar43 + ai22*ai33*ai44 + + - ai22*ai34*ai43 - ai22*ar11*ar33 - ai22*ar11*ar44 - ai22*ar33*ar44 + + + ai22*ar34*ar43 - ai24*ai33*ai42 + ai24*ar11*ar42 - ai24*ar32*ar43 + + + ai24*ar33*ar42 - ai33*ar11*ar22 - ai33*ar11*ar44 - ai33*ar22*ar44 + + + ai33*ar24*ar42 + ai34*ar11*ar43 + ai34*ar22*ar43 - ai34*ar23*ar42 + + + ai42*ar11*ar24 - ai42*ar23*ar34 + ai42*ar24*ar33 + ai43*ar11*ar34 + + + ai43*ar22*ar34 - ai43*ar24*ar32 - ai44*ar11*ar22 - ai44*ar11*ar33 + + - ai44*ar22*ar33 + ai44*ar23*ar32 + +cpr02 := - ai11*ai22 - ai11*ai33 - ai11*ai44 - ai22*ai33 - ai22*ai44 + + + ai24*ai42 - ai33*ai44 + ai34*ai43 + ar11*ar22 + ar11*ar33 + + + ar11*ar44 + ar22*ar33 + ar22*ar44 - ar23*ar32 - ar24*ar42 + + + ar33*ar44 - ar34*ar43 + +cpi02 := ai11*ar22 + ai11*ar33 + ai11*ar44 + ai22*ar11 + ai22*ar33 + ai22*ar44 + + - ai24*ar42 + ai33*ar11 + ai33*ar22 + ai33*ar44 - ai34*ar43 + + - ai42*ar24 - ai43*ar34 + ai44*ar11 + ai44*ar22 + ai44*ar33 + +cpr03 := - (ar11 + ar22 + ar33 + ar44) + +cpi03 := - (ai11 + ai22 + ai33 + ai44) + + 2 2 +rpr00 := cpi00 + cpr00 + +rpr01 := 2*(cpi00*cpi01 + cpr00*cpr01) + + 2 2 +rpr02 := 2*cpi00*cpi02 + cpi01 + 2*cpr00*cpr02 + cpr01 + +rpr03 := 2*(cpi00*cpi03 + cpi01*cpi02 + cpr00*cpr03 + cpr01*cpr02) + + 2 2 +rpr04 := 2*cpi01*cpi03 + cpi02 + 2*cpr00 + 2*cpr01*cpr03 + cpr02 + +rpr05 := 2*(cpi02*cpi03 + cpr01 + cpr02*cpr03) + + 2 2 +rpr06 := cpi03 + 2*cpr02 + cpr03 + +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 - (ai*i + ar) + +x2:=lam-(br+i*bi); + + +x2 := lam - (bi*i + br) + +x3:=lam-(cr+i*ci); + + +x3 := lam - (ci*i + cr) + +hurwitzp x1; + +Necessary and sufficient conditions are: + - ar > 0 + + +cond + + +% Example H.2 + +x:=hurw(x0*x1); + + +x := 2*lam*( - ai*i - ar + 1) + 2*(ai*i + ar + 1) + +hurwitzp x; + +Necessary and sufficient conditions are: + 2 2 +4*( - ai - ar + 1) > 0 + + +cond + + +% Example H.3 + +x:=(x1*x2); + + + 2 +x := lam - lam*(ai*i + ar + bi*i + br) - ai*bi + ai*br*i + ar*bi*i + ar*br + +hurwitzp x; + +Necessary and sufficient conditions are: + - (ar + br) > 0 + + 2 2 2 2 +ar*br*(ai - 2*ai*bi + ar + 2*ar*br + bi + br ) > 0 + + +cond + + +% Example H.4 + +x:=(x1*x2*x3); + + + 3 2 +x := lam - lam *(ai*i + ar + bi*i + br + ci*i + cr) + lam*( - ai*bi + ai*br*i + + - ai*ci + ai*cr*i + ar*bi*i + ar*br + ar*ci*i + ar*cr - bi*ci + bi*cr*i + + + br*ci*i + br*cr) + ai*bi*ci*i + ai*bi*cr + ai*br*ci - ai*br*cr*i + + + ar*bi*ci - ar*bi*cr*i - ar*br*ci*i - ar*br*cr + +hurwitzp x; + +Necessary and sufficient conditions are: + - (ar + br + cr) > 0 + + 2 2 3 3 +ai *ar*br + ai *ar*cr - 2*ai*ar*bi*br - 2*ai*ar*ci*cr + ar *br + ar *cr + + 2 2 2 2 2 2 3 2 + + 2*ar *br + 4*ar *br*cr + 2*ar *cr + ar*bi *br + ar*br + 4*ar*br *cr + + 2 2 3 2 3 + + 4*ar*br*cr + ar*ci *cr + ar*cr + bi *br*cr - 2*bi*br*ci*cr + br *cr + + 2 2 2 3 + + 2*br *cr + br*ci *cr + br*cr > 0 + + 4 2 4 4 2 4 4 2 4 2 +ar*br*cr*( - ai *bi + 2*ai *bi*ci - ai *br - 2*ai *br*cr - ai *ci - ai *cr + + 3 3 3 2 3 2 3 + + 2*ai *bi - 2*ai *bi *ci + 2*ai *bi*br + 4*ai *bi*br*cr + + 3 2 3 2 3 2 3 + - 2*ai *bi*ci + 2*ai *bi*cr + 2*ai *br *ci + 4*ai *br*ci*cr + + 3 3 3 2 2 2 2 2 2 + + 2*ai *ci + 2*ai *ci*cr - 2*ai *ar *bi + 4*ai *ar *bi*ci + + 2 2 2 2 2 2 2 2 2 2 2 + - 2*ai *ar *br - 4*ai *ar *br*cr - 2*ai *ar *ci - 2*ai *ar *cr + + 2 2 2 2 2 + - 2*ai *ar*bi *br - 2*ai *ar*bi *cr + 4*ai *ar*bi*br*ci + + 2 2 3 2 2 + + 4*ai *ar*bi*ci*cr - 2*ai *ar*br - 6*ai *ar*br *cr + + 2 2 2 2 2 2 2 3 + - 2*ai *ar*br*ci - 6*ai *ar*br*cr - 2*ai *ar*ci *cr - 2*ai *ar*cr + + 2 4 2 3 2 2 2 2 2 + - ai *bi - 2*ai *bi *ci - 2*ai *bi *br - 2*ai *bi *br*cr + + 2 2 2 2 2 2 2 2 2 + + 6*ai *bi *ci - 2*ai *bi *cr - 2*ai *bi*br *ci - 8*ai *bi*br*ci*cr + + 2 3 2 2 2 4 2 3 + - 2*ai *bi*ci - 2*ai *bi*ci*cr - ai *br - 2*ai *br *cr + + 2 2 2 2 2 2 2 2 2 3 + - 2*ai *br *ci - 2*ai *br *cr - 2*ai *br*ci *cr - 2*ai *br*cr + + 2 4 2 2 2 2 4 2 3 2 2 + - ai *ci - 2*ai *ci *cr - ai *cr + 2*ai*ar *bi - 2*ai*ar *bi *ci + + 2 2 2 2 2 + + 2*ai*ar *bi*br + 4*ai*ar *bi*br*cr - 2*ai*ar *bi*ci + + 2 2 2 2 2 + + 2*ai*ar *bi*cr + 2*ai*ar *br *ci + 4*ai*ar *br*ci*cr + + 2 3 2 2 3 2 + + 2*ai*ar *ci + 2*ai*ar *ci*cr + 4*ai*ar*bi *cr + 4*ai*ar*bi *br*ci + + 2 2 2 + - 8*ai*ar*bi *ci*cr + 4*ai*ar*bi*br *cr - 8*ai*ar*bi*br*ci + + 2 2 3 + + 8*ai*ar*bi*br*cr + 4*ai*ar*bi*ci *cr + 4*ai*ar*bi*cr + + 3 2 3 + + 4*ai*ar*br *ci + 8*ai*ar*br *ci*cr + 4*ai*ar*br*ci + + 2 4 3 2 3 2 + + 4*ai*ar*br*ci*cr + 2*ai*bi *ci - 2*ai*bi *ci + 2*ai*bi *cr + + 2 2 2 2 3 + + 4*ai*bi *br *ci + 4*ai*bi *br*ci*cr - 2*ai*bi *ci + + 2 2 2 2 2 2 + - 2*ai*bi *ci*cr - 2*ai*bi*br *ci + 2*ai*bi*br *cr + + 2 3 4 2 2 + + 4*ai*bi*br*ci *cr + 4*ai*bi*br*cr + 2*ai*bi*ci + 4*ai*bi*ci *cr + + 4 4 3 2 3 + + 2*ai*bi*cr + 2*ai*br *ci + 4*ai*br *ci*cr + 2*ai*br *ci + + 2 2 4 2 4 4 2 4 + + 2*ai*br *ci*cr - ar *bi + 2*ar *bi*ci - ar *br - 2*ar *br*cr + + 4 2 4 2 3 2 3 2 3 + - ar *ci - ar *cr - 2*ar *bi *br - 2*ar *bi *cr + 4*ar *bi*br*ci + + 3 3 3 3 2 3 2 + + 4*ar *bi*ci*cr - 2*ar *br - 6*ar *br *cr - 2*ar *br*ci + + 3 2 3 2 3 3 2 4 2 3 + - 6*ar *br*cr - 2*ar *ci *cr - 2*ar *cr - ar *bi + 2*ar *bi *ci + + 2 2 2 2 2 2 2 2 2 2 2 + - 2*ar *bi *br - 6*ar *bi *br*cr - 2*ar *bi *ci - 2*ar *bi *cr + + 2 2 2 2 3 + + 2*ar *bi*br *ci + 8*ar *bi*br*ci*cr + 2*ar *bi*ci + + 2 2 2 4 2 3 2 2 2 + + 2*ar *bi*ci*cr - ar *br - 6*ar *br *cr - 2*ar *br *ci + + 2 2 2 2 2 2 3 2 4 + - 10*ar *br *cr - 6*ar *br*ci *cr - 6*ar *br*cr - ar *ci + + 2 2 2 2 4 4 3 + - 2*ar *ci *cr - ar *cr - 2*ar*bi *cr + 4*ar*bi *ci*cr + + 2 2 2 2 2 2 + - 4*ar*bi *br *cr - 2*ar*bi *br*ci - 6*ar*bi *br*cr + + 2 2 2 3 2 3 + - 2*ar*bi *ci *cr - 2*ar*bi *cr + 4*ar*bi*br *ci*cr + 4*ar*bi*br*ci + + 2 4 3 2 3 2 + + 4*ar*bi*br*ci*cr - 2*ar*br *cr - 2*ar*br *ci - 6*ar*br *cr + + 2 2 2 3 4 2 2 + - 6*ar*br *ci *cr - 6*ar*br *cr - 2*ar*br*ci - 4*ar*br*ci *cr + + 4 4 2 4 2 3 3 3 2 + - 2*ar*br*cr - bi *ci - bi *cr + 2*bi *ci + 2*bi *ci*cr + + 2 2 2 2 2 2 2 2 2 3 + - 2*bi *br *ci - 2*bi *br *cr - 2*bi *br*ci *cr - 2*bi *br*cr + + 2 4 2 2 2 2 4 2 3 2 2 + - bi *ci - 2*bi *ci *cr - bi *cr + 2*bi*br *ci + 2*bi*br *ci*cr + + 4 2 4 2 3 2 3 3 2 4 + - br *ci - br *cr - 2*br *ci *cr - 2*br *cr - br *ci + + 2 2 2 2 4 + - 2*br *ci *cr - br *cr ) > 0 + + +cond + +clear x,x0,x1,x2,x3; + + + +%*********************************************************************** +%***** ***** +%***** T e s t Examples --- Module L I N B A N D ***** +%***** ***** +%*********************************************************************** + +on evallhseqp; + + % So both sides of equations evaluate. + +% 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=0.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.0e+000 + 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.0e+000 + 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.0 + arer=0.0 + 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=0.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.0e+000 + 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.0e+000 + 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.0 + arer=0.0 + 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=0.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.0e+000 + acof(n,2)=1.0 + acof(n,3)=0.0e+000 + 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.0e+000 + acof(n,2)=1.0 + acof(n,3)=0.0e+000 + 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.0 + arer=0.0 + 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=0.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.0e+000 + 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.0e+000 + 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.0 + arer=0.0 + 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 v(x..one),w(x..half); + + +iim aa, v, diff(v,t)=c*diff(w,x), + w, diff(w,t)=c*diff(v,x); + + +***************************** +***** Program ***** IIMET Ver 1.1.2 +***************************** + + 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 one grid point + Equation for w variable is integrated in half grid point + + Equations after Discretization Using IIM : + ========================================== + + 2*j - 1 2*j + 1 2*j + 1 2*j - 1 +((w(n,---------) - w(n,---------) - w(n + 1,---------) + w(n + 1,---------))*c + 2 2 2 2 + + *ht + 2*(v(n + 1,j) - v(n,j))*hx)/(2*ht*hx) = 0 + + +( - (v(n,j + 1) - v(n,j) - v(n + 1,j) + v(n + 1,j + 1))*c*ht + + 2*j + 1 2*j + 1 + + 2*(w(n + 1,---------) - w(n,---------))*hx)/(2*ht*hx) = 0 + 2 2 + + + +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 + + 2 +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 + + 2 +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 +a(1,1) := (2*v(n + 1,j)*hx - 2*v(n,j)*hx + w(n + 1,---------)*c*ht + 2 + + 2*j + 1 2*j - 1 + - w(n + 1,---------)*c*ht + w(n,---------)*c*ht + 2 2 + + 2*j + 1 + - w(n,---------)*c*ht)/(2*ht*hx) + 2 + +a(1,2):=aa(1,0); + + +a(1,2) := ( - v(n + 1,j + 1)*c*ht + v(n + 1,j)*c*ht - v(n,j + 1)*c*ht + + 2*j + 1 2*j + 1 + + v(n,j)*c*ht + 2*w(n + 1,---------)*hx - 2*w(n,---------)*hx)/(2*ht + 2 2 + + *hx) + +off prfourmat; + + +b:=ampmat a; + + + kx*hx +ax := ------- + 2 + + + [ 2 2 2 2 ] + [ - sin(ax) *c *ht + hx 2*i*sin(ax)*c*ht*hx ] + [-------------------------- ----------------------- ] + [ 2 2 2 2 2 2 2 2 ] + [ sin(ax) *c *ht + hx sin(ax) *c *ht + hx ] +b := [ ] + [ 2 2 2 2 ] + [ 2*i*sin(ax)*c*ht*hx - 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 *ht *hx + 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 2 2 2 + + 2*sin(ax) *c *ht *hx + hx )/(sin(ax) *c *ht + 2*sin(ax) *c *ht *hx + + 4 + + 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 *ht *hx + 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 2 2 2 + + 2*sin(ax) *c *ht *hx + hx )/(sin(ax) *c *ht + 2*sin(ax) *c *ht *hx + + 4 + + hx ) + +pol:=hurw num pol; + + +pol := + + 2 2 2 2 2 2 2 2 2 2 2 2 2 +4*lam *sin(ax) *c *ht *(sin(ax) *c *ht + hx ) + 4*hx *(sin(ax) *c *ht + hx ) + +hurwitzp pol; + +Necessary and sufficient conditions are: + 2 + hx +----------------- > 0 + 2 2 2 + sin(ax) *c *ht + + +cond + +bt:=tcon b; + + + [ 2 2 2 2 ] + [ - sin(ax) *c *ht + hx - 2*sin(ax)*c*ht*hx*i ] + [-------------------------- ------------------------ ] + [ 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*ht*hx*i - 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*i*sin(ax)*c*ht ] + [ 1 ------------------ ] + [ hx ] + [ ] +b := [ 2 2 2 2 ] + [ 2*i*sin(ax)*c*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 + + +cond + +bt:=tcon b; + + + [ - 2*sin(ax)*c*ht*i ] + [ 1 --------------------- ] + [ hx ] + [ ] +bt := [ 2 2 2 2 ] + [ - 2*sin(ax)*c*ht*i - 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 *ht *i ] +[------------------------- --------------------- ] +[ 2 3 ] +[ hx hx ] +[ ] +[ 3 3 3 4 4 4 2 2 2 2 4 ] +[ - 8*sin(ax) *c *ht *i 16*sin(ax) *c *ht - 4*sin(ax) *c *ht *hx + hx ] +[------------------------ --------------------------------------------------] +[ 3 4 ] +[ hx hx ] + + +bt*b-b*bt; + + +[ 3 3 3 ] +[ 16*sin(ax) *c *ht *i ] +[ 0 ----------------------] +[ 3 ] +[ hx ] +[ ] +[ 3 3 3 ] +[ - 16*sin(ax) *c *ht *i ] +[------------------------- 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 - i*sin(ax)*c*ht +b := mat((-------------------------------------------,-------------------, + hx hx + + - i*sin(ay)*c*ht + -------------------), + hx + + - i*sin(ax)*c*ht cos(ax)*c*ht - c*ht + hx + (-------------------,--------------------------,0), + hx hx + + - i*sin(ay)*c*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) + + 2 2 2 2 + + lam*hx*(3*cos(ax)*cos(ay)*c *ht - 5*cos(ax)*c *ht + + 2 2 + + 4*cos(ax)*c*ht*hx - 5*cos(ay)*c *ht + 4*cos(ay)*c*ht*hx + + 2 2 2 3 3 + + 7*c *ht - 8*c*ht*hx + 3*hx ) + 4*cos(ax)*cos(ay)*c *ht + + 2 2 3 3 2 2 + - 3*cos(ax)*cos(ay)*c *ht *hx - 4*cos(ax)*c *ht + 5*cos(ax)*c *ht *hx + + 2 3 3 2 2 + - 2*cos(ax)*c*ht*hx - 4*cos(ay)*c *ht + 5*cos(ay)*c *ht *hx + + 2 3 3 2 2 2 3 3 + - 2*cos(ay)*c*ht*hx + 4*c *ht - 7*c *ht *hx + 4*c*ht*hx - 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 +pol := lam + lam *(4*s1 *cap + 4*s2 *cap - 3) + lam + + 2 2 2 2 2 2 2 2 2 + *(12*s1 *s2 *cap + 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*lam*cap + + 2 2 2 2 2 2 2 2 2 + *(3*s1 *s2 *cap - 3*s1 *s2 *cap - s1 *cap + s1 - s2 *cap + s2 ) + 8*( + + 2 2 3 2 2 2 2 2 2 2 2 + - 2*s1 *s2 *cap + 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 + + 2 2 2 2 2 2 2 2 2 + (3) 16*cap*(3*s1 *s2 *cap - 3*s1 *s2 *cap - s1 *cap + s1 - s2 *cap + s2 ) + + > 0 + + 2 2 3 2 2 2 2 2 2 2 2 + (4) 8*( - 2*s1 *s2 *cap + 3*s1 *s2 *cap + s1 *cap - 2*s1 *cap + s2 *cap + + 2 + - 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 4 4 + + 8*s1 *s2 *cap - 10*s1 *s2 *cap + 3*s1 *s2 - s1 *cap + s1 + + 2 4 2 2 4 2 4 2 2 + + 8*s1 *s2 *cap - 10*s1 *s2 *cap + 3*s1 *s2 - 2*s1 *s2 *cap + + 2 2 4 4 + + s1 *s2 - s2 *cap + 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 6 4 4 + - 27*s1 *s2 *cap - 32*s1 *s2 *cap + 100*s1 *s2 *cap + + 6 4 3 6 4 2 6 2 4 + - 93*s1 *s2 *cap + 27*s1 *s2 *cap + 10*s1 *s2 *cap + + 6 2 3 6 2 2 6 2 6 3 + - 31*s1 *s2 *cap + 26*s1 *s2 *cap - 6*s1 *s2 *cap - s1 *cap + + 6 2 6 4 6 5 4 6 4 + + 3*s1 *cap - 2*s1 *cap - 32*s1 *s2 *cap + 100*s1 *s2 *cap + + 4 6 3 4 6 2 4 4 4 + - 93*s1 *s2 *cap + 27*s1 *s2 *cap + 20*s1 *s2 *cap + + 4 4 3 4 4 2 4 4 + - 76*s1 *s2 *cap + 73*s1 *s2 *cap - 21*s1 *s2 *cap + + 4 2 3 4 2 2 4 2 + - 3*s1 *s2 *cap + 16*s1 *s2 *cap - 14*s1 *s2 *cap + + 4 2 4 4 2 6 4 + + 3*s1 *s2 - s1 *cap + s1 + 10*s1 *s2 *cap + + 2 6 3 2 6 2 2 6 + - 31*s1 *s2 *cap + 26*s1 *s2 *cap - 6*s1 *s2 *cap + + 2 4 3 2 4 2 2 4 + - 3*s1 *s2 *cap + 16*s1 *s2 *cap - 14*s1 *s2 *cap + + 2 4 2 2 2 2 6 3 6 2 + + 3*s1 *s2 - 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 sin(ax)*c*ht*i +bt := mat((-------------------------------------------,----------------, + hx hx + + sin(ay)*c*ht*i + ----------------), + hx + + sin(ax)*c*ht*i cos(ax)*c*ht - c*ht + hx + (----------------,--------------------------,0), + hx hx + + sin(ay)*c*ht*i 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*ht*hx + + 2 2 2 2 2 2 + - 4*cos(ay)*c *ht + 2*cos(ay)*c*ht*hx + 6*c *ht - 4*c*ht*hx + hx )/hx , + + 2 2 2 2 + sin(ax)*c *ht *i*( - cos(ay) + 1) sin(ay)*c *ht *i*( - cos(ax) + 1) + -----------------------------------,-----------------------------------), + 2 2 + hx hx + + 2 2 + sin(ax)*c *ht *i*(cos(ay) - 1) + (--------------------------------, + 2 + hx + + 2 2 2 2 2 + - 2*cos(ax)*c *ht + 2*cos(ax)*c*ht*hx + 2*c *ht - 2*c*ht*hx + hx + ----------------------------------------------------------------------, + 2 + hx + + 2 2 + sin(ax)*sin(ay)*c *ht + ------------------------), + 2 + hx + + 2 2 2 2 + sin(ay)*c *ht *i*(cos(ax) - 1) sin(ax)*sin(ay)*c *ht + (--------------------------------,------------------------, + 2 2 + hx hx + + 2 2 2 2 2 + - 2*cos(ay)*c *ht + 2*cos(ay)*c*ht*hx + 2*c *ht - 2*c*ht*hx + hx + ----------------------------------------------------------------------)) + 2 + hx + + +bt*b-b*bt; + + + 2 2 + 2*sin(ax)*c *ht *i*( - cos(ay) + 1) +mat((0,-------------------------------------, + 2 + hx + + 2 2 + 2*sin(ay)*c *ht *i*( - cos(ax) + 1) + -------------------------------------), + 2 + hx + + 2 2 + 2*sin(ax)*c *ht *i*(cos(ay) - 1) + (----------------------------------,0,0), + 2 + hx + + 2 2 + 2*sin(ay)*c *ht *i*(cos(ax) - 1) + (----------------------------------,0,0)) + 2 + hx + + +clear a,b,bt,pol; + + + +%*********************************************************************** + +end; +(TIME: fide 15030 15900)