Codemist Standard Lisp 3.54 for DEC Alpha: May 23 1994
Dump file created: Mon May 23 10:39:11 1994
REDUCE 3.5, 15-Oct-93 ...
Memory allocation: 6023424 bytes
+++ About to read file tstlib.red
sf(1) := 1
sf(2) := 1
sf(3) := 1
%***********************************************************************
%***** *****
%***** Package F I D E - Test Examples Ver. 1.1.1 November 2,1993 *****
%***** *****
%***********************************************************************
%***********************************************************************
%***** *****
%***** T e s t Examples --- Module E X P R E S *****
%***** *****
%***********************************************************************
let cos th**2=1 - sin th**2,
cos fi**2=1 - sin fi**2;
factor df;
on rat;
for all x,y let diff(x,y)=df(x,y);
depend u,r,th,fi;
depend v,r,th,fi;
depend f,r,th,fi;
depend w,r,th,fi;
% Spherical coordinate system
scalefactors 3,r*sin th*cos fi,r*sin th*sin fi,r*cos th,r,th,fi;
sf(1) := 1
sf(2) := r
sf(3) := sin(th)*r
tensor a1,a2,a3,a4,a5;
vectors u,v;
dyads w;
a1:=grad f;
a1:= ( df(f,r) ,
df(f,th)
---------- ,
r
df(f,fi)
----------- )
sin(th)*r
a2:=div u;
df(u(3),fi) df(u(2),th)
a2:=------------- + ------------- + df(u(1),r)
sin(th)*r r
cos(th)*u(2) + 2*sin(th)*u(1)
+ -------------------------------
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:= ( ------------------ + ------------------ + ---------------
2 2 2 2
sin(th)*r r sin(th) *r
2*df(v(1),r) df(v(1),th,2)
+ df(v(1),r,2) + -------------- + ---------------
r 2
r
df(v(1),th)*cos(th) 2*(cos(th)*v(2) + sin(th)*v(1))
+ --------------------- - ---------------------------------
2 2
sin(th)*r 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),r,2) + --------------
2 2 r
sin(th) *r
df(v(3),th,2) df(v(3),th)*cos(th)
+ --------------- + ---------------------
2 2
r sin(th)*r
2*df(v(2),fi)*cos(th) 2*df(v(1),fi) - v(3)
+ ----------------------- + --------------- + -------------
2 2 2 2 2
sin(th) *r sin(th)*r sin(th) *r
)
a3:=2*a3+a4;
- 2*df(v(3),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)
+ ------------------ + --------------- + df(v(1),r,2)
2 2 2
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(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),r,2) + --------------
2 2 r
sin(th) *r
df(v(3),th,2) df(v(3),th)*cos(th)
+ --------------- + ---------------------
2 2
r sin(th)*r
2*df(v(2),fi)*cos(th) 2*df(v(1),fi)
+ ----------------------- + 2*df(v(2),r) + ---------------
2 2 2
sin(th) *r sin(th)*r
2
- 2*df(v(1),th) 2*sin(th) *v(2)*r - v(3)
+ ------------------ + -------------------------- )
r 2 2
sin(th) *r
a5:=lapl f;
df(f,fi,2) 2*df(f,r) df(f,th,2)
a5:=------------- + df(f,r,2) + ----------- + ------------
2 2 r 2
sin(th) *r r
df(f,th)*cos(th)
+ ------------------
2
sin(th)*r
a1:=a1+div w;
df(w(3,1),fi) df(w(2,1),th)
a1:= ( --------------- + --------------- + df(w(1,1),r) + df(f,r) + (
sin(th)*r r
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;
sf(1) := 1
sf(2) := 1
sf(3) := 1
%***********************************************************************
%***** *****
%***** T e s t Examples --- Module I I M E T *****
%***** *****
%***********************************************************************
% Example I.1 - 1-D Lagrangian Hydrodynamics
off exp;
factor diff;
on rat,eqfu;
% Declare which indexes will be given to coordinates
coordinates x,t into j,m;
% Declares uniform grid in x coordinate
grid uniform,x;
% Declares dependencies of functions on coordinates
dependence eta(t,x),v(t,x),eps(t,x),p(t,x);
% Declares p as known function
given p;
same eta,v,p;
iim a, eta,diff(eta,t)-eta*diff(v,x)=0,
v,diff(v,t)+eta/ro*diff(p,x)=0,
eps,diff(eps,t)+eta*p/ro*diff(v,x)=0;
*****************************
***** Program ***** IIMET Ver 1.1 17/05/91
*****************************
Partial Differential Equations
==============================
diff(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 :
==========================================
( - ((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*(eta(j + 1,m + 1) - eta(j + 1,m) + eta(j,m + 1) - eta(j,m))*hx)
/(4*(ht(m + 1) + ht(m))*hx) = 0
(((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*(v(j + 1,m + 1) - v(j + 1,m) + v(j,m + 1) - v(j,m))*hx*ro)/(4
*(ht(m + 1) + ht(m))*hx*ro) = 0
(((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
*(eps(j + 1,m + 1) - eps(j + 1,m) + eps(j,m + 1) - eps(j,m))*hx*ro)/
(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 17/05/91
*****************************
Partial Differential Equations
==============================
diff(u,t) + diff(u,x) + diff(v,x,2)*cos(x) = 0
diff(sin(x)*u,x) + diff(v,t) = 0
0 interpolations are needed in x coordinate
Equation for u variable is integrated in half grid point
Equation for v variable is integrated in half grid point
0 interpolations are needed in t coordinate
Equation for u variable is integrated in half grid point
Equation for v variable is integrated in half grid point
Equations after Discretization Using IIM :
==========================================
((u(j + 1,m + 1) + u(j + 1,m) - u(j,m + 1) - u(j,m))*ht*hx
2
+ (u(j + 1,m + 1) - u(j + 1,m) + u(j,m + 1) - u(j,m))*hx + (
2*j - 1 2*j - 1 2*j + 3
v(---------,m + 1) + v(---------,m) + v(---------,m + 1)
2 2 2
2*j + 3 2*j + 1 2*j + 1
+ v(---------,m) - 2*v(---------,m + 1) - 2*v(---------,m))
2 2 2
2*j + 1 2
*cos(x(---------))*ht)/(2*ht*hx ) = 0
2
((u(j + 1,m + 1) + u(j + 1,m))*sin(x(j + 1))*ht
- (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*ht*hx) = 0
2 2
clear a;
%***********************************************************************
% Example I.3 - Schrodinger equation
factor diff;
coordinates t,x into m,j;
grid uniform,x,t;
dependence ur(x,t),ui(x,t);
same ui,ur;
iim a,ur,-diff(ui,t)+1/2*diff(ur,x,2)+(ur**2+ui**2)*ur=0,
ui,diff(ur,t)+1/2*diff(ui,x,2)+(ur**2+ui**2)*ui=0;
*****************************
***** Program ***** IIMET Ver 1.1 17/05/91
*****************************
Partial Differential Equations
==============================
diff(ur,x,2) 2 2
- diff(ui,t) + -------------- = - ur*(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 + 1,j - 1) + ur(m + 1,j + 1) - 2*ur(m + 1,j) + ur(m,j - 1)
2 2
+ ur(m,j + 1) - 2*ur(m,j))*ht - 4*ui(m + 1,j)*hx + 4*ui(m,j)*hx )
2 2 2
/(4*ht*hx ) = - (ui(m + 1,j) *ur(m + 1,j) + ui(m,j) *ur(m,j)
3 3
+ ur(m + 1,j) + ur(m,j) )/2
((ui(m + 1,j - 1) + ui(m + 1,j + 1) - 2*ui(m + 1,j) + ui(m,j - 1)
2 2
+ ui(m,j + 1) - 2*ui(m,j))*ht + 4*ur(m + 1,j)*hx - 4*ur(m,j)*hx )
2 3 2
/(4*ht*hx ) = - (ui(m + 1,j) + ui(m + 1,j)*ur(m + 1,j)
3 2
+ ui(m,j) + ui(m,j)*ur(m,j) )/2
clear a;
clearsame;
%***********************************************************************
% Example I.4 - Vector calculus in p.d.e. input
% cooperation with expres module
% 2-D hydrodynamics
scalefactors 2,x,y,x,y;
sf(1) := 1
sf(2) := 1
vectors u;
off exp,twogrid;
on eqfu;
factor diff,ht,hx,hy;
coordinates x,y,t into j,i,m;
grid uniform,x,y,t;
dependence n(t,x,y),u(t,x,y),p(t,x,y);
iim a,n,diff(n,t)+u.grad n+n*div u=0,
u,m*n*(diff(u,t)+u.grad u)+grad p=vect(0,0),
p,3/2*(diff(p,t)+u.grad p)+5/2*p*div u=0;
*****************************
***** Program ***** IIMET Ver 1.1 17/05/91
*****************************
Partial Differential Equations
==============================
diff(n,t) + diff(n,x)*u1 + diff(n,y)*u2 + diff(u1,x)*n + diff(u2,y)*n
= 0
diff(p,x) + diff(u1,t)*m*n + diff(u1,x)*m*n*u1 + diff(u1,y)*m*n*u2
= 0
diff(p,y) + diff(u2,t)*m*n + diff(u2,x)*m*n*u1 + diff(u2,y)*m*n*u2
= 0
3*diff(p,t) 3*diff(p,x)*u1 3*diff(p,y)*u2
------------- + ---------------- + ----------------
2 2 2
5*diff(u1,x)*p 5*diff(u2,y)*p
+ ---------------- + ---------------- = 0
2 2
0 interpolations are needed in x coordinate
Equation for n variable is integrated in half grid point
Equation for u1 variable is integrated in half grid point
Equation for u2 variable is integrated in half grid point
Equation for p variable is integrated in half grid point
0 interpolations are needed in y coordinate
Equation for n variable is integrated in half grid point
Equation for u1 variable is integrated in half grid point
Equation for u2 variable is integrated in half grid point
Equation for p variable is integrated in half grid point
0 interpolations are needed in t coordinate
Equation for n variable is integrated in half grid point
Equation for u1 variable is integrated in half grid point
Equation for u2 variable is integrated in half grid point
Equation for p variable is integrated in half grid point
Equations after Discretization Using IIM :
==========================================
-1 -1
(hy *hx *((u1(j + 1,i + 1,m + 1)*hy + u2(j + 1,i + 1,m + 1)*hx)
*n(j + 1,i + 1,m + 1) +
(u1(j + 1,i + 1,m)*hy + u2(j + 1,i + 1,m)*hx)
*n(j + 1,i + 1,m) +
(u1(j + 1,i,m + 1)*hy - u2(j + 1,i,m + 1)*hx)
*n(j + 1,i,m + 1)
+ (u1(j + 1,i,m)*hy - u2(j + 1,i,m)*hx)*n(j + 1,i,m) -
(u1(j,i + 1,m + 1)*hy - u2(j,i + 1,m + 1)*hx)
*n(j,i + 1,m + 1)
- (u1(j,i + 1,m)*hy - u2(j,i + 1,m)*hx)*n(j,i + 1,m)
- (u1(j,i,m + 1)*hy + u2(j,i,m + 1)*hx)*n(j,i,m + 1)
-1
- (u1(j,i,m)*hy + u2(j,i,m)*hx)*n(j,i,m)))/4 + (ht *(
n(j + 1,i + 1,m + 1) - n(j + 1,i + 1,m) + n(j + 1,i,m + 1)
- n(j + 1,i,m) + n(j,i + 1,m + 1) - n(j,i + 1,m)
+ n(j,i,m + 1) - n(j,i,m)))/4 = 0
-1
(hx *((p(j + 1,i + 1,m + 1) + p(j + 1,i + 1,m) + p(j + 1,i,m + 1)
+ p(j + 1,i,m))
- (p(j,i + 1,m + 1) + p(j,i + 1,m) + p(j,i,m + 1) + p(j,i,m))
-1 -1
))/4 + (hy *hx *(((
n(j + 1,i + 1,m + 1)*u1(j + 1,i + 1,m + 1)
+ n(j,i + 1,m + 1)*u1(j,i + 1,m + 1))
*(u1(j + 1,i + 1,m + 1) - u1(j,i + 1,m + 1)) + (
n(j + 1,i + 1,m)*u1(j + 1,i + 1,m)
+ n(j,i + 1,m)*u1(j,i + 1,m))
*(u1(j + 1,i + 1,m) - u1(j,i + 1,m)) + (
n(j + 1,i,m + 1)*u1(j + 1,i,m + 1)
+ n(j,i,m + 1)*u1(j,i,m + 1))
*(u1(j + 1,i,m + 1) - u1(j,i,m + 1)) + (
n(j + 1,i,m)*u1(j + 1,i,m)
+ n(j,i,m)*u1(j,i,m))
*(u1(j + 1,i,m) - u1(j,i,m)))*hy + ((
n(j + 1,i + 1,m + 1)*u2(j + 1,i + 1,m + 1)
+ n(j + 1,i,m + 1)*u2(j + 1,i,m + 1))
*(u1(j + 1,i + 1,m + 1) - u1(j + 1,i,m + 1)) + (
n(j + 1,i + 1,m)*u2(j + 1,i + 1,m)
+ n(j + 1,i,m)*u2(j + 1,i,m))
*(u1(j + 1,i + 1,m) - u1(j + 1,i,m)) + (
n(j,i + 1,m + 1)*u2(j,i + 1,m + 1)
+ n(j,i,m + 1)*u2(j,i,m + 1))
*(u1(j,i + 1,m + 1) - u1(j,i,m + 1)) + (
n(j,i + 1,m)*u2(j,i + 1,m)
+ n(j,i,m)*u2(j,i,m))
-1
*(u1(j,i + 1,m) - u1(j,i,m)))*hx)*m)/8 + (ht *(
(n(j + 1,i + 1,m + 1) + n(j + 1,i + 1,m))
*(u1(j + 1,i + 1,m + 1) - u1(j + 1,i + 1,m)) +
(n(j + 1,i,m + 1) + n(j + 1,i,m))
*(u1(j + 1,i,m + 1) - u1(j + 1,i,m)) +
(n(j,i + 1,m + 1) + n(j,i + 1,m))
*(u1(j,i + 1,m + 1) - u1(j,i + 1,m))
+ (n(j,i,m + 1) + n(j,i,m))*(u1(j,i,m + 1) - u1(j,i,m)))*m)/8
= 0
-1
(hy *((n(j + 1,i + 1,m + 1)*u2(j + 1,i + 1,m + 1)
+ n(j + 1,i,m + 1)*u2(j + 1,i,m + 1))
*(u2(j + 1,i + 1,m + 1) - u2(j + 1,i,m + 1)) + (
n(j + 1,i + 1,m)*u2(j + 1,i + 1,m)
+ n(j + 1,i,m)*u2(j + 1,i,m))
*(u2(j + 1,i + 1,m) - u2(j + 1,i,m)) + (
n(j,i + 1,m + 1)*u2(j,i + 1,m + 1)
+ n(j,i,m + 1)*u2(j,i,m + 1))
*(u2(j,i + 1,m + 1) - u2(j,i,m + 1)) +
(n(j,i + 1,m)*u2(j,i + 1,m) + n(j,i,m)*u2(j,i,m))
-1
*(u2(j,i + 1,m) - u2(j,i,m)))*m)/8 + (hx *((
n(j + 1,i + 1,m + 1)*u1(j + 1,i + 1,m + 1)
+ n(j,i + 1,m + 1)*u1(j,i + 1,m + 1))
*(u2(j + 1,i + 1,m + 1) - u2(j,i + 1,m + 1)) + (
n(j + 1,i + 1,m)*u1(j + 1,i + 1,m)
+ n(j,i + 1,m)*u1(j,i + 1,m))
*(u2(j + 1,i + 1,m) - u2(j,i + 1,m)) + (
n(j + 1,i,m + 1)*u1(j + 1,i,m + 1)
+ n(j,i,m + 1)*u1(j,i,m + 1))
*(u2(j + 1,i,m + 1) - u2(j,i,m + 1)) +
(n(j + 1,i,m)*u1(j + 1,i,m) + n(j,i,m)*u1(j,i,m))
-1 -1
*(u2(j + 1,i,m) - u2(j,i,m)))*m)/8 + (hy *ht *((
(n(j + 1,i + 1,m + 1) + n(j + 1,i + 1,m))
*(u2(j + 1,i + 1,m + 1) - u2(j + 1,i + 1,m)) +
(n(j + 1,i,m + 1) + n(j + 1,i,m))
*(u2(j + 1,i,m + 1) - u2(j + 1,i,m)) +
(n(j,i + 1,m + 1) + n(j,i + 1,m))
*(u2(j,i + 1,m + 1) - u2(j,i + 1,m))
+ (n(j,i,m + 1) + n(j,i,m))*(u2(j,i,m + 1) - u2(j,i,m)))*hy
*m + 2*(p(j + 1,i + 1,m + 1) + p(j + 1,i + 1,m))*ht
- 2*(p(j + 1,i,m + 1) + p(j + 1,i,m))*ht
+ 2*(p(j,i + 1,m + 1) + p(j,i + 1,m))*ht
- 2*(p(j,i,m + 1) + p(j,i,m))*ht))/8 = 0
-1 -1
(hy *hx *(5*((p(j + 1,i + 1,m + 1) + p(j + 1,i,m + 1))
*(u2(j + 1,i + 1,m + 1) - u2(j + 1,i,m + 1)) +
(p(j + 1,i + 1,m) + p(j + 1,i,m))
*(u2(j + 1,i + 1,m) - u2(j + 1,i,m)) +
(p(j,i + 1,m + 1) + p(j,i,m + 1))
*(u2(j,i + 1,m + 1) - u2(j,i,m + 1)) +
(p(j,i + 1,m) + p(j,i,m))*(u2(j,i + 1,m) - u2(j,i,m)))
*hx + 3*((p(j + 1,i + 1,m + 1) - p(j + 1,i,m + 1))
*(u2(j + 1,i + 1,m + 1) + u2(j + 1,i,m + 1)) +
(p(j + 1,i + 1,m) - p(j + 1,i,m))
*(u2(j + 1,i + 1,m) + u2(j + 1,i,m)) +
(p(j,i + 1,m + 1) - p(j,i,m + 1))
*(u2(j,i + 1,m + 1) + u2(j,i,m + 1)) +
(p(j,i + 1,m) - p(j,i,m))*(u2(j,i + 1,m) + u2(j,i,m)))
*hx + 5*((p(j + 1,i + 1,m + 1) + p(j,i + 1,m + 1))
*(u1(j + 1,i + 1,m + 1) - u1(j,i + 1,m + 1)) +
(p(j + 1,i + 1,m) + p(j,i + 1,m))
*(u1(j + 1,i + 1,m) - u1(j,i + 1,m)) +
(p(j + 1,i,m + 1) + p(j,i,m + 1))
*(u1(j + 1,i,m + 1) - u1(j,i,m + 1)) +
(p(j + 1,i,m) + p(j,i,m))*(u1(j + 1,i,m) - u1(j,i,m)))
*hy + 3*((p(j + 1,i + 1,m + 1) - p(j,i + 1,m + 1))
*(u1(j + 1,i + 1,m + 1) + u1(j,i + 1,m + 1)) +
(p(j + 1,i + 1,m) - p(j,i + 1,m))
*(u1(j + 1,i + 1,m) + u1(j,i + 1,m)) +
(p(j + 1,i,m + 1) - p(j,i,m + 1))
*(u1(j + 1,i,m + 1) + u1(j,i,m + 1)) +
(p(j + 1,i,m) - p(j,i,m))*(u1(j + 1,i,m) + u1(j,i,m)))
-1
*hy))/16 + (3*ht *(p(j + 1,i + 1,m + 1)
- p(j + 1,i + 1,m) + p(j + 1,i,m + 1) - p(j + 1,i,m)
+ p(j,i + 1,m + 1) - p(j,i + 1,m) + p(j,i,m + 1) - p(j,i,m)))/
8 = 0
clear a,u;
%***********************************************************************
% Example I.5 - 1-D hydrodynamics up to 3-rd moments (heat flow)
coordinates x,t into j,m;
grid uniform,x,t;
dependence n(x,t),u(x,t),tt(x,t),p(x,t),q(x,t);
iim a, n,diff(n,t)+u*diff(n,x)+diff(u,x)=0,
u,n*m*(diff(u,t)+u*diff(u,x))+k*diff(n*tt,x)+diff(p,x)=0,
tt,3/2*k*n*(diff(tt,t)+u*diff(tt,x))+n*k*tt*diff(u,x)+1/2*p
*diff(u,x)+diff(q,x)=0,
p,diff(p,t)+u*diff(p,x)+p*diff(u,x)+n*k*tt*diff(u,x)+2/5*diff(q,x)
=0,
q,diff(q,t)+u*diff(q,x)+q*diff(u,x)+5/2*n*k**2*tt/m*diff(tt,x)+n*k
*tt*diff(p,x)-p*diff(p,x)=0;
*****************************
***** Program ***** IIMET Ver 1.1 17/05/91
*****************************
Partial Differential Equations
==============================
diff(n,t) + diff(n,x)*u + diff(u,x) = 0
diff(n*tt,x)*k + diff(p,x) + diff(u,t)*m*n + diff(u,x)*m*n*u =
0
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)
5
= 0
diff(p,x)*(k*n*tt - p) + diff(q,t) + diff(q,x)*u
2
5*diff(tt,x)*k *n*tt
+ ---------------------- + diff(u,x)*q = 0
2*m
0 interpolations are needed in x coordinate
Equation for n variable is integrated in half grid point
Equation for u variable is integrated in half grid point
Equation for tt variable is integrated in half grid point
Equation for p variable is integrated in half grid point
Equation for q variable is integrated in half grid point
0 interpolations are needed in t coordinate
Equation for n variable is integrated in half grid point
Equation for u variable is integrated in half grid point
Equation for tt variable is integrated in half grid point
Equation for p variable is integrated in half grid point
Equation for q variable is integrated in half grid point
Equations after Discretization Using IIM :
==========================================
-1
(hx *(((n(j + 1,m + 1) - n(j,m + 1))*(u(j + 1,m + 1) + u(j,m + 1))
+ (n(j + 1,m) - n(j,m))*(u(j + 1,m) + u(j,m)))
+ 2*(u(j + 1,m + 1) + u(j + 1,m)) - 2*(u(j,m + 1) + u(j,m))))
-1
ht *(n(j + 1,m + 1) - n(j + 1,m) + n(j,m + 1) - n(j,m))
/4 + ----------------------------------------------------------
2
= 0
-1
(hx *(((n(j + 1,m + 1)*u(j + 1,m + 1) + n(j,m + 1)*u(j,m + 1))
*(u(j + 1,m + 1) - u(j,m + 1)) +
(n(j + 1,m)*u(j + 1,m) + n(j,m)*u(j,m))*(u(j + 1,m) - u(j,m))
)*m + 2*(n(j + 1,m + 1)*tt(j + 1,m + 1)
+ n(j + 1,m)*tt(j + 1,m) - n(j,m + 1)*tt(j,m + 1)
- n(j,m)*tt(j,m))*k
+ 2*(p(j + 1,m + 1) + p(j + 1,m)) - 2*(p(j,m + 1) + p(j,m))))
-1
/4 + (ht *((n(j + 1,m + 1) + n(j + 1,m))
*(u(j + 1,m + 1) - u(j + 1,m))
+ (n(j,m + 1) + n(j,m))*(u(j,m + 1) - u(j,m)))*m)/4
= 0
-1
(hx *(2*((n(j + 1,m + 1)*tt(j + 1,m + 1) + n(j,m + 1)*tt(j,m + 1))
*(u(j + 1,m + 1) - u(j,m + 1)) +
(n(j + 1,m)*tt(j + 1,m) + n(j,m)*tt(j,m))
*(u(j + 1,m) - u(j,m)))*k + 3*(
(n(j + 1,m + 1)*u(j + 1,m + 1) + n(j,m + 1)*u(j,m + 1))
*(tt(j + 1,m + 1) - tt(j,m + 1)) +
(n(j + 1,m)*u(j + 1,m) + n(j,m)*u(j,m))
*(tt(j + 1,m) - tt(j,m)))*k + (
(p(j + 1,m + 1) + p(j,m + 1))*(u(j + 1,m + 1) - u(j,m + 1))
+ (p(j + 1,m) + p(j,m))*(u(j + 1,m) - u(j,m)))
+ 4*(q(j + 1,m + 1) + q(j + 1,m)) - 4*(q(j,m + 1) + q(j,m))))
-1
/8 + (3*ht *((n(j + 1,m + 1) + n(j + 1,m))
*(tt(j + 1,m + 1) - tt(j + 1,m))
+ (n(j,m + 1) + n(j,m))*(tt(j,m + 1) - tt(j,m)))*k)/8
= 0
-1
(hx *(5*((n(j + 1,m + 1)*tt(j + 1,m + 1) + n(j,m + 1)*tt(j,m + 1))
*(u(j + 1,m + 1) - u(j,m + 1)) +
(n(j + 1,m)*tt(j + 1,m) + n(j,m)*tt(j,m))
*(u(j + 1,m) - u(j,m)))*k + 5*(
(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))) + 5*(
(p(j + 1,m + 1) - p(j,m + 1))*(u(j + 1,m + 1) + u(j,m + 1))
+ (p(j + 1,m) - p(j,m))*(u(j + 1,m) + u(j,m)))
+ 4*(q(j + 1,m + 1) + q(j + 1,m)) - 4*(q(j,m + 1) + q(j,m))))
-1
ht *(p(j + 1,m + 1) - p(j + 1,m) + p(j,m + 1) - p(j,m))
/20 + ----------------------------------------------------------
2
= 0
-1
(hx *(2*((n(j + 1,m + 1)*tt(j + 1,m + 1) + n(j,m + 1)*tt(j,m + 1))
*(p(j + 1,m + 1) - p(j,m + 1)) +
(n(j + 1,m)*tt(j + 1,m) + n(j,m)*tt(j,m))
*(p(j + 1,m) - p(j,m)))*k*m + 5*(
(n(j + 1,m + 1)*tt(j + 1,m + 1) + n(j,m + 1)*tt(j,m + 1))
*(tt(j + 1,m + 1) - tt(j,m + 1)) +
(n(j + 1,m)*tt(j + 1,m) + n(j,m)*tt(j,m))
2
*(tt(j + 1,m) - tt(j,m)))*k - 2*(
(p(j + 1,m + 1) + p(j,m + 1))*(p(j + 1,m + 1) - p(j,m + 1))
+ (p(j + 1,m) + p(j,m))*(p(j + 1,m) - p(j,m)))*m + 2*(
(q(j + 1,m + 1) + q(j,m + 1))*(u(j + 1,m + 1) - u(j,m + 1))
+ (q(j + 1,m) + q(j,m))*(u(j + 1,m) - u(j,m)))*m + 2*(
(q(j + 1,m + 1) - q(j,m + 1))*(u(j + 1,m + 1) + u(j,m + 1))
+ (q(j + 1,m) - q(j,m))*(u(j + 1,m) + u(j,m)))*m))/(8*m)
-1
ht *(q(j + 1,m + 1) - q(j + 1,m) + q(j,m + 1) - q(j,m))
+ ---------------------------------------------------------- = 0
2
clear a;
remfac diff,ht,hx,hy;
on exp;
off rat;
%***********************************************************************
%***** *****
%***** T e s t Examples --- Module A P P R O X *****
%***** *****
%***********************************************************************
% Example A.1
coordinates x,t into j,n;
maxorder t=2,x=3;
functions u,v;
approx( (u(n+1/2)-u(n-1/2))/ht=(v(n+1/2,j+1/2)-v(n+1/2,j-1/2)
+v(n-1/2,j+1/2)-v(n-1/2,j-1/2))/(2*hx) );
Difference scheme approximates differential equation df(u,t)=df(v,x)
with orders of approximation:
2
hx
2
ht
% Example A.2
maxorder t=3,x=3;
approx( (u(n+1)-u(n))/ht=(u(n+1,j+1/2)-u(n+1,j-1/2)
+u(n,j+1/2)-u(n,j-1/2))/(2*hx) );
Difference scheme approximates differential equation df(u,t)=df(u,x)
with orders of approximation:
2
hx
ht
% Example A.3
maxorder t=2,x=3;
center t=1/2;
approx( (u(n+1)-u(n))/ht=(v(n+1,j+1/2)-v(n+1,j-1/2)
+v(n,j+1/2)-v(n,j-1/2))/(2*hx) );
Difference scheme approximates differential equation df(u,t)=df(v,x)
with orders of approximation:
2
hx
2
ht
% Example A.4
approx( u(n+1)/ht=(v(n+1,j+1/2)-v(n+1,j-1/2)
+v(n,j+1/2)-v(n,j-1/2))/(2*hx) );
Reformulate difference scheme, grid steps remain in denominators
Difference scheme approximates differential equation 0=df(v,x)
with orders of approximation:
2
hx
-1
ht
% Example A.5
maxorder t=3,x=3;
approx( (u(n+1)-u(n))/ht=(u(n+1,j+1/2)-u(n+1,j-1/2))/hx);
Difference scheme approximates differential equation df(u,t)=df(u,x)
with orders of approximation:
2
hx
ht
% Example A.6
approx( (u(n+1)-u(n))/ht=(u(n+1/2,j+1/2)-u(n+1/2,j-1/2))/hx);
Difference scheme approximates differential equation df(u,t)=df(u,x)
with orders of approximation:
2
hx
2
ht
% Example A.7;
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)
((sin(ax)*ht
2 2
*(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
ai21 := (sin(ax)*ht*(cos(ax)*sin(ax) *c0 *ht *u0
2 2 2 2 2
- cos(ax)*c0 *ht*hx + 2*sin(ax) *c0 *ht *hx*u0
2 2 2 3
+ 2*sin(ax) *ht *hx*r0*u0 - 2*hx *r0))/(2*(
4 4 4 3 3 3
sin(ax) *ht *u0 - 2*sin(ax) *ht *hx*i*u0
2 2 2 2 3 4
- 2*sin(ax) *ht *hx *u0 + 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
ai22 := (sin(ax)*ht*( - cos(ax)*sin(ax) *c0 *ht *u0
2 2 2 2 2
+ cos(ax)*c0 *ht*hx + 2*sin(ax) *c0 *ht *hx*u0
2 2 3 3
- 2*sin(ax) *ht *hx*u0 - 2*hx *u0))/(2*(
4 4 4 3 3 3
sin(ax) *ht *u0 - 2*sin(ax) *ht *hx*i*u0
2 2 2 2 3 4
- 2*sin(ax) *ht *hx *u0 + 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
bb := mat(((cos(ax)*c*ht*hy + cos(ay)*c*ht*hx - c*ht*hx - c*ht*hy
- i*sin(ax)*c*ht - i*sin(ay)*c*ht
+ hx*hy)/(hx*hy),-------------------,-------------------
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) +
2 2 2 2
lam*(3*cos(ax)*cos(ay)*c *ht *hx*hy - 3*cos(ax)*c *ht *hx*hy
2 2 2 2
- 2*cos(ax)*c *ht *hy + 4*cos(ax)*c*ht*hx*hy
2 2 2 2 2 2 2
+ cos(ay) *c *ht *hx - 2*cos(ay)*c *ht *hx
2 2 2
- 3*cos(ay)*c *ht *hx*hy + 4*cos(ay)*c*ht*hx *hy
2 2 2 2 2 2 2 2 2
+ sin(ay) *c *ht *hx + c *ht *hx + 3*c *ht *hx*hy
2 2 2 2 2
+ 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
- 3*cos(ax)*cos(ay)*c *ht *hx*hy
2 3 3 3 3
- cos(ax)*sin(ay) *c *ht *hx - cos(ax)*c *ht *hx
3 3 2 2
- 2*cos(ax)*c *ht *hy + 3*cos(ax)*c *ht *hx*hy
2 2 2 2
+ 2*cos(ax)*c *ht *hy - 2*cos(ax)*c*ht*hx*hy
2 3 3 2 2 2 2
+ cos(ay) *c *ht *hx - cos(ay) *c *ht *hx
3 3 3 3
- 2*cos(ay)*c *ht *hx - 2*cos(ay)*c *ht *hy
2 2 2 2 2
+ 2*cos(ay)*c *ht *hx + 3*cos(ay)*c *ht *hx*hy
2 2 3 3
- 2*cos(ay)*c*ht*hx *hy + sin(ay) *c *ht *hx
2 2 2 2 3 3 3 3 2 2 2
- sin(ay) *c *ht *hx + c *ht *hx + 2*c *ht *hy - c *ht *hx
2 2 2 2 2 2
- 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
pol := lam + lam *(4*s1 *cap1 + 4*s2 *cap2 - 3) + lam*(
2 2 2 2 2
12*s1 *s2 *cap1*cap2 + 4*s1 *cap1 - 8*s1 *cap1
2 2 2 2 2 2
+ 4*s2 *cap2 - 8*s2 *cap2 + 3) + 8*s1 *s2 *cap1 *cap2
2 2 2 2 2 2 2
+ 8*s1 *s2 *cap1*cap2 - 12*s1 *s2 *cap1*cap2 - 4*s1 *cap1
2 2 2 2
+ 4*s1 *cap1 - 4*s2 *cap2 + 4*s2 *cap2 - 1
clear cos ax,cos ay,sin ax,sin ay,cos ax2**2,cos ay2**2,sin ax2,sin ay2,
hx,hy;
pol:=complexpol pol;
2 2
If 8*s1 *s2 *cap1*cap2*(cap1 + cap2) = 0 and 0
= 0 , a root of the polynomial is equal to 1.
3 2 2 2
pol := lam + lam *(4*s1 *cap1 + 4*s2 *cap2 - 3) + lam*(
2 2 2 2 2
12*s1 *s2 *cap1*cap2 + 4*s1 *cap1 - 8*s1 *cap1
2 2 2 2 2 2
+ 4*s2 *cap2 - 8*s2 *cap2 + 3) + 8*s1 *s2 *cap1 *cap2
2 2 2 2 2 2 2
+ 8*s1 *s2 *cap1*cap2 - 12*s1 *s2 *cap1*cap2 - 4*s1 *cap1
2 2 2 2
+ 4*s1 *cap1 - 4*s2 *cap2 + 4*s2 *cap2 - 1
pol1:=hurw pol;
3 2 2 2
pol1 := 8*lam *s1 *s2 *cap1*cap2*(cap1 + cap2) + 8*lam *(
2 2 2 2 2 2
- 3*s1 *s2 *cap1 *cap2 - 3*s1 *s2 *cap1*cap2
2 2 2 2 2 2
+ 3*s1 *s2 *cap1*cap2 + s1 *cap1 + s2 *cap2 ) + 8*lam*(
2 2 2 2 2 2
3*s1 *s2 *cap1 *cap2 + 3*s1 *s2 *cap1*cap2
2 2 2 2 2
- 6*s1 *s2 *cap1*cap2 - 2*s1 *cap1 + 2*s1 *cap1
2 2 2 2 2 2
- 2*s2 *cap2 + 2*s2 *cap2) + 8*( - s1 *s2 *cap1 *cap2
2 2 2 2 2 2 2
- s1 *s2 *cap1*cap2 + 3*s1 *s2 *cap1*cap2 + s1 *cap1
2 2 2 2
- 2*s1 *cap1 + s2 *cap2 - 2*s2 *cap2 + 1)
denotid cp;
pol:=denotepol pol;
3 2
pol := lam + lam *cpr02 + lam*cpr01 + cpr00
prdenot;
2 2 2 2 2 2
cpr00 := 8*s1 *s2 *cap1 *cap2 + 8*s1 *s2 *cap1*cap2
2 2 2 2 2
- 12*s1 *s2 *cap1*cap2 - 4*s1 *cap1 + 4*s1 *cap1
2 2 2
- 4*s2 *cap2 + 4*s2 *cap2 - 1
2 2 2 2 2
cpr01 := 12*s1 *s2 *cap1*cap2 + 4*s1 *cap1 - 8*s1 *cap1
2 2 2
+ 4*s2 *cap2 - 8*s2 *cap2 + 3
2 2
cpr02 := 4*s1 *cap1 + 4*s2 *cap2 - 3
cleardenot;
clear aa,bb,pol,pol1;
%***********************************************************************
% Example C.6 : Lax-Wendrov (V. Ganzha)
coordinates t,x,y into n,m,k;
grid uniform,t,x,y;
let cos ax**2=1-sin ax**2,
cos ay**2=1-sin ay**2;
unfunc u1,u2,u3,u4;
matrix aa(1,4),bb(4,4);
aa(1,1):=4*(u1(n+1)-u1(n))/ht+
(w*(u1(m+2)-u1(m-2)+u1(m+1,k+1)+u1(m+1,k-1)-
u1(m-1,k+1)-u1(m-1,k-1))+p*(u2(m+2)-u2(m-2)+u2(m+1,k+1)+
u2(m+1,k-1)-u2(m-1,k+1)-u2(m-1,k-1))+
v*(u1(m+1,k+1)+u1(m-1,k+1)-
u1(m+1,k-1)-u1(m-1,k-1)+u1(k+2)-u1(k-2))+p*(u3(m+1,k+1)+
u3(m-1,k+1)-u3(m+1,k-1)-u3(m-1,k-1)+u3(k+2)-u3(k-2)))/hx+ht*
(2*w**2*(-u1(m+2)+2*u1(m)-u1(m-2))+4*w*p*(-u2(m+2)+2*u2(m)-
u2(m-2))+2*(-u4(m+2)+2*u4(m)-u4(m-2))+2*v**2*(-u1(k+2)+
2*u1(k)-u1(k-2))+4*v*p*(u3(k+2)+2*u3(k)-u3(k-2))+2*(-u4(k+2)+
2*u4(k)-u4(k-2))+4*w*v*(-u1(m+1,k+1)+u1(m+1,k-1)+u1(m-1,k+1)-
u1(m-1,k-1))+4*p*v*(-u2(m+1,k+1)+u2(m+1,k-1)+u2(m-1,k+1)-
u2(m-1,k-1))+4*w*p*(-u3(m+1,k+1)+u3(m+1,k-1)+u3(m-1,k+1)-
u3(m-1,k-1)))/hx/hx$
aa(1,2):=4*p*(u2(n+1)-u2(n))/ht+
(w*p*(u2(m+2)-u2(m-2)+u2(m+1,k+1)+
u2(m+1,k-1)-u2(m-1,k+1)-u2(m-1,k-1))+u4(m+2)-u4(m-2)+
u4(m+1,k+1)+
u4(m+1,k-1)-u4(m-1,k+1)-u4(m-1,k-1)+
p*v*(u2(m+1,k+1)+u2(m-1,k+1)+
u2(k+2)-u2(k-2)-u2(m+1,k-1)-u2(m-1,k-1)))/hx+ht*(2*w**2*p*
(-u2(m+2)+2*u2(m)-u2(m-2))+2*p*c**2*(-u2(m+2)+2*u2(m)-u2(m-2))
+4*w*(-u4(m+2)+2*u4(m)-u4(m-2))+2*p*v**2*(-u2(k+2)+2*u2(k)-
u2(k-2))+4*w*p*v*(-u2(m+1,k+1)+u2(m+1,k-1)+u2(m-1,k+1)-
u2(m-1,k-1))+2*p*c**2*(-u3(m+1,k+1)+u3(m+1,k-1)+u3(m-1,k+1)
-u3(m-1,k-1))+4*v*(-u4(m+1,k+1)+u4(m+1,k-1)+u4(m-1,k+1)-
u4(m-1,k-1)))/hx/hx$
aa(1,3):=4*p*(u3(n+1)-u3(n))/ht+(w*p*(u3(m+2)-u3(m-2)+u3(m+1,k+1)+
u3(m+1,k-1)-u3(m-1,k+1)-u3(m-1,k-1))+u4(k+2)-u4(k-2)+
u4(m+1,k+1)-u4(m+1,k-1)+u4(m-1,k+1)-u4(m-1,k-1)+
p*v*(u3(m+1,k+1)+u3(m-1,k+1)+u3(k+2)-u3(k-2)-u3(m+1,k-1)-
u3(m-1,k-1)))/hx+ht*(2*w**2*p*(-u3(m+2)+2*u3(m)-u3(m-2))+
2*p*c**2*(-u3(k+2)+2*u3(k)-u3(k-2))+4*v*(-u4(k+2)+
2*u4(k)-u4(k-2))+2*p*v**2*(-u3(k+2)+2*u3(k)-u3(k-2))+
4*w*p*v*(-u3(m+1,k+1)+u3(m+1,k-1)+u3(m-1,k+1)-
u3(m-1,k-1))+2*p*c**2*(-u2(m+1,k+1)+u2(m+1,k-1)+
u2(m-1,k+1)-u2(m-1,k-1))+4*w*(u4(m+1,k+1)+u4(m+1,k-1)+
u4(m-1,k+1)-u4(m-1,k-1)))/hx/hx$
aa(1,4):=4*(u4(n+1)-u4(n))/ht+(p*c**2*(u2(m+2)-u2(m-2)+u2(m+1,k+1)+
u2(m+1,k-1)-u2(m-1,k+1)-u2(m-1,k-1))+w*(u4(m+2)-
u4(m-2)+u4(m+1,k+1)+u4(m+1,k-1)-u4(m-1,k+1)-u4(m-1,k-1))+
+p*c**2*(u3(m+1,k+1)+u3(m-1,k+1)-u3(m+1,k-1)-
u3(m-1,k-1)+u3(k+2)-u3(k-2))+v*(u4(m+1,k+1)+u4(m-1,k+1)-
u4(m+1,k-1)-u4(m-1,k-1)+u4(k+2)-u4(k-2)))/hx+ht*
(2*w**2*(-u4(m+2)+2*u4(m)-u4(m-2))+4*w*p*c**2*(-u2(m+2)+
2*u2(m)-u2(m-2))+2*c**2*(-u4(m+2)+2*u4(m)-u4(m-2))+
4*p*v*c**2*(-u3(k+2)+2*u3(k)-u3(k-2))+2*c**2*(-u4(k+2)+
2*u4(k)-u4(k-2))+2*v**2*(-u4(k+2)+2*u4(k)-u4(k-2))+
4*p*v*c**2*(-u2(m+1,k+1)+u2(m+1,k-1)+u2(m-1,k+1)-
u2(m-1,k-1))+4*w*p*c**2*(-u3(m+1,k+1)+u3(m+1,k-1)+
u3(m-1,k+1)-u3(m-1,k-1))+4*w*v*(-u4(m+1,k+1)+
u4(m+1,k-1)+u4(m-1,k+1)-u4(m-1,k-1)))/hx/hx$
bb:=ampmat aa;
ax := kx*hx
ay := ky*hy
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*sin(ax) *ht *w - 4*sin(ax)*sin(ay)*ht *v*w
2 2 2 2 2
- 2*sin(ay) *ht *v + hx )/hx ,(sin(ax)*ht*p*(
- i*cos(ax)*hx - i*cos(ay)*hx - 4*sin(ax)*ht*w
2
- 4*sin(ay)*ht*v))/hx ,(ht*p*(
- i*cos(ax)*sin(ay)*hx - 4*i*cos(ay)*sin(ay)*ht*v
- i*cos(ay)*sin(ay)*hx - 4*sin(ax)*sin(ay)*ht*w
2 2 2
2 - 2*ht *(sin(ax) + sin(ay) )
- 2*ht*v))/hx ,--------------------------------),
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
2 2 2
- i*cos(ay)*sin(ay)*ht*hx*v - 2*sin(ax) *c *ht
2 2 2 2
- 2*sin(ax) *ht *w - 4*sin(ax)*sin(ay)*ht *v*w
2 2 2 2 2
- 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
- i*cos(ay)*hx - 4*sin(ax)*ht*w - 4*sin(ay)*ht*v))/
2
(hx *p)),
2 2
- 2*sin(ax)*sin(ay)*c *ht
(0,-----------------------------,(
2
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
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*sin(ax)*sin(ay)*ht*w
2 2
- 4*sin(ay) *ht*v))/(hx *p)),
2
(0,(sin(ax)*c *ht*p*( - i*cos(ax)*hx - i*cos(ay)*hx
2 2
- 4*sin(ax)*ht*w - 4*sin(ay)*ht*v))/hx ,(sin(ay)*c
*ht*p*( - i*cos(ax)*hx - i*cos(ay)*hx - 4*sin(ax)*ht*w
2
- 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
2 2 2
- i*cos(ay)*sin(ay)*ht*hx*v - 2*sin(ax) *c *ht
2 2 2 2
- 2*sin(ax) *ht *w - 4*sin(ax)*sin(ay)*ht *v*w
2 2 2 2 2 2 2 2
- 2*sin(ay) *c *ht - 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
2
- ar22 - ar33 - ar44) + 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 4
+ 2*lam *(cpi02*cpi03 + cpr01 + cpr02*cpr03) + lam
2 2
*(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
pol := lam + lam *rpr07 + lam *rpr06 + lam *rpr05 + lam *rpr04
3 2
+ lam *rpr03 + lam *rpr02 + lam*rpr01 + rpr00
prdenot;
2 2 2 2
ar11 := - 2*s1 *k1 - 4*s1*s2*k1*k2 - 2*s2 *k2 + 1
ai11 := - (s1*c1*k1 + s1*c2*k1 + s2*c1*k2 + s2*c2*k2)
ar12 := - 4*s1*p*r1*(s1*k1 + s2*k2)
ai12 := - s1*p*r1*(c1 + c2)
ar13 := 2*p*r1*( - 2*s1*s2*k1 - k2)
ai13 := s2*p*r1*( - c1 - 4*c2*k2 - c2)
2 2 2
ar14 := - 2*r1 *(s1 + s2 )
2 2 2 2 2 2
ar22 := - 2*s1 *k1 - 2*s1 *k3 - 4*s1*s2*k1*k2 - 2*s2 *k2 + 1
ai22 := - (s1*c1*k1 + s1*c2*k1 + s2*c1*k2 + s2*c2*k2)
2
ar23 := - 2*s1*s2*k3
- 4*s1*r1*(s1*k1 + s2*k2)
ar24 := ----------------------------
p
- s1*r1*(c1 + c2)
ai24 := --------------------
p
2
ar32 := - 2*s1*s2*k3
2 2 2 2 2 2
ar33 := - 2*s1 *k1 - 4*s1*s2*k1*k2 - 2*s2 *k2 - 2*s2 *k3 + 1
ai33 := - (s1*c1*k1 + s1*c2*k1 + s2*c1*k2 + s2*c2*k2)
2
2*r1*( - s1*s2*k1 - 2*s2 *k2 - c1*c2*k1)
ar34 := ------------------------------------------
p
r1*( - 2*s1*c2*k1 - 2*s2*c1*k1 - s2*c1 - s2*c2)
ai34 := -------------------------------------------------
p
2
- 4*s1*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
ar44 :=
2 2 2 2 2 2 2 2
- 2*s1 *k1 - 2*s1 *k3 - 4*s1*s2*k1*k2 - 2*s2 *k2 - 2*s2 *k3 + 1
ai44 := - (s1*c1*k1 + s1*c2*k1 + s2*c1*k2 + s2*c2*k2)
cpr00 := 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
ai *ar*br + ai *ar*cr - 2*ai*ar*bi*br - 2*ai*ar*ci*cr + ar *br
3 2 2 2 2 2 2 3
+ ar *cr + 2*ar *br + 4*ar *br*cr + 2*ar *cr + ar*bi *br + ar*br
2 2 2 3 2
+ 4*ar*br *cr + 4*ar*br*cr + ar*ci *cr + ar*cr + bi *br*cr
3 2 2 2 3
- 2*bi*br*ci*cr + br *cr + 2*br *cr + br*ci *cr + br*cr > 0
4 2 4 4 2 4 4 2
ar*br*cr*( - ai *bi + 2*ai *bi*ci - ai *br - 2*ai *br*cr - ai *ci
4 2 3 3 3 2 3 2
- ai *cr + 2*ai *bi - 2*ai *bi *ci + 2*ai *bi*br
3 3 2 3 2
+ 4*ai *bi*br*cr - 2*ai *bi*ci + 2*ai *bi*cr
3 2 3 3 3 3 2
+ 2*ai *br *ci + 4*ai *br*ci*cr + 2*ai *ci + 2*ai *ci*cr
2 2 2 2 2 2 2 2
- 2*ai *ar *bi + 4*ai *ar *bi*ci - 2*ai *ar *br
2 2 2 2 2 2 2 2
- 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*ai *ar*br*ci - 6*ai *ar*br*cr - 2*ai *ar*ci *cr
2 3 2 4 2 3 2 2 2
- 2*ai *ar*cr - ai *bi - 2*ai *bi *ci - 2*ai *bi *br
2 2 2 2 2 2 2 2
- 2*ai *bi *br*cr + 6*ai *bi *ci - 2*ai *bi *cr
2 2 2 2 3
- 2*ai *bi*br *ci - 8*ai *bi*br*ci*cr - 2*ai *bi*ci
2 2 2 4 2 3 2 2 2
- 2*ai *bi*ci*cr - ai *br - 2*ai *br *cr - 2*ai *br *ci
2 2 2 2 2 2 3 2 4
- 2*ai *br *cr - 2*ai *br*ci *cr - 2*ai *br*cr - ai *ci
2 2 2 2 4 2 3 2 2
- 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*ai*ar *ci + 2*ai*ar *ci*cr + 4*ai*ar*bi *cr
2 2
+ 4*ai*ar*bi *br*ci - 8*ai*ar*bi *ci*cr
2 2
+ 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
+ 4*ai*ar*br*ci*cr + 2*ai*bi *ci - 2*ai*bi *ci
3 2 2 2 2
+ 2*ai*bi *cr + 4*ai*bi *br *ci + 4*ai*bi *br*ci*cr
2 3 2 2 2 2
- 2*ai*bi *ci - 2*ai*bi *ci*cr - 2*ai*bi*br *ci
2 2 2 3
+ 2*ai*bi*br *cr + 4*ai*bi*br*ci *cr + 4*ai*bi*br*cr
4 2 2 4
+ 2*ai*bi*ci + 4*ai*bi*ci *cr + 2*ai*bi*cr
4 3 2 3
+ 2*ai*br *ci + 4*ai*br *ci*cr + 2*ai*br *ci
2 2 4 2 4 4 2
+ 2*ai*br *ci*cr - ar *bi + 2*ar *bi*ci - ar *br
4 4 2 4 2 3 2
- 2*ar *br*cr - ar *ci - ar *cr - 2*ar *bi *br
3 2 3 3
- 2*ar *bi *cr + 4*ar *bi*br*ci + 4*ar *bi*ci*cr
3 3 3 2 3 2 3 2
- 2*ar *br - 6*ar *br *cr - 2*ar *br*ci - 6*ar *br*cr
3 2 3 3 2 4 2 3
- 2*ar *ci *cr - 2*ar *cr - ar *bi + 2*ar *bi *ci
2 2 2 2 2 2 2 2
- 2*ar *bi *br - 6*ar *bi *br*cr - 2*ar *bi *ci
2 2 2 2 2 2
- 2*ar *bi *cr + 2*ar *bi*br *ci + 8*ar *bi*br*ci*cr
2 3 2 2 2 4 2 3
+ 2*ar *bi*ci + 2*ar *bi*ci*cr - ar *br - 6*ar *br *cr
2 2 2 2 2 2 2 2
- 2*ar *br *ci - 10*ar *br *cr - 6*ar *br*ci *cr
2 3 2 4 2 2 2 2 4
- 6*ar *br*cr - ar *ci - 2*ar *ci *cr - ar *cr
4 3 2 2
- 2*ar*bi *cr + 4*ar*bi *ci*cr - 4*ar*bi *br *cr
2 2 2 2 2 2
- 2*ar*bi *br*ci - 6*ar*bi *br*cr - 2*ar*bi *ci *cr
2 3 2 3
- 2*ar*bi *cr + 4*ar*bi*br *ci*cr + 4*ar*bi*br*ci
2 4 3 2
+ 4*ar*bi*br*ci*cr - 2*ar*br *cr - 2*ar*br *ci
3 2 2 2 2 3
- 6*ar*br *cr - 6*ar*br *ci *cr - 6*ar*br *cr
4 2 2 4 4 2
- 2*ar*br*ci - 4*ar*br*ci *cr - 2*ar*br*cr - bi *ci
4 2 3 3 3 2 2 2 2
- bi *cr + 2*bi *ci + 2*bi *ci*cr - 2*bi *br *ci
2 2 2 2 2 2 3 2 4
- 2*bi *br *cr - 2*bi *br*ci *cr - 2*bi *br*cr - bi *ci
2 2 2 2 4 2 3 2 2
- 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.0
arhs(n)=v(1)
n=n+1
alcd(n)=1.0
ad(n)=-2.0
aucd(n)=1.0
arhs(n)=v(3)-(2.0*v(2))+v(1)
DO 25002 k=3,99,1
n=n+1
alcd(n)=1.0
ad(n)=-2.0
aucd(n)=1.0
arhs(n)=v(k-1)+v(k+1)-(2.0*v(k))
25002 CONTINUE
n=n+1
alcd(n)=1.0
ad(n)=-2.0
aucd(n)=1.0
arhs(n)=v(101)-(2.0*v(100))+v(99)
n=n+1
alcd(n)=0.0
ad(n)=1.0
arhs(n)=v(101)
CALL DGTSL(N,ALCD,AD,AUCD,ARHS,IER)
C N IS NUMBER OF EQUATIONS
C ALCD,AD,AUCD,ARHS ARE ARRAYS OF DIMENSION AT LEAST (N)
C IF (IER.NE.0) MATRIX ACOF IS ALGORITHMICALLY SINGULAR
IF(IER.NE.0) CALL ERROUT
n=1
u(1)=arhs(n)
n=n+1
u(2)=arhs(n)
DO 25003 k=3,99,1
n=n+1
u(k)=arhs(n)
25003 CONTINUE
n=n+1
u(100)=arhs(n)
n=n+1
u(101)=arhs(n)
amer=0.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.0
arhs(n)=v(1)
n=n+1
alcd(n)=1.0
ad(n)=-2.0
aucd(n+1)=1.0
arhs(n)=v(3)-(2.0*v(2))+v(1)
DO 25006 k=3,99,1
n=n+1
alcd(n)=1.0
ad(n)=-2.0
aucd(n+1)=1.0
arhs(n)=v(k-1)+v(k+1)-(2.0*v(k))
25006 CONTINUE
n=n+1
alcd(n)=1.0
ad(n)=-2.0
aucd(n+1)=1.0
arhs(n)=v(101)-(2.0*v(100))+v(99)
n=n+1
alcd(n)=0.0
ad(n)=1.0
arhs(n)=v(101)
ier=0
CALL F01LEF(N,AD,0.,AUCD,ALCD,1.E-10,AU2CD,IN,IER)
C N IS NUMBER OF EQUATIONS
C ALCD,AD,AUCD,AU2CD,ARHS ARE ARRAYS OF DIMENSION AT LEAST (N)
C IN IS INTEGER ARRAY OF DIMENSION AT LEAST (N)
IF(IER.NE.0 .OR. IN(N).NE.0) CALL ERROUT
CALL F04LEF(1,N,AD,AUCD,ALCD,AU2CD,IN,ARHS,0.,IER)
IF(IER.NE.0) CALL ERROUT
n=1
u(1)=arhs(n)
n=n+1
u(2)=arhs(n)
DO 25007 k=3,99,1
n=n+1
u(k)=arhs(n)
25007 CONTINUE
n=n+1
u(100)=arhs(n)
n=n+1
u(101)=arhs(n)
amer=0.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.0
acof(n,2)=1.0
acof(n,3)=0.0
arhs(n)=v(1)
n=n+1
acof(n,1)=1.0
acof(n,2)=-2.0
acof(n,3)=1.0
arhs(n)=v(3)-(2.0*v(2))+v(1)
DO 25010 k=3,99,1
n=n+1
acof(n,1)=1.0
acof(n,2)=-2.0
acof(n,3)=1.0
arhs(n)=v(k-1)+v(k+1)-(2.0*v(k))
25010 CONTINUE
n=n+1
acof(n,1)=1.0
acof(n,2)=-2.0
acof(n,3)=1.0
arhs(n)=v(101)-(2.0*v(100))+v(99)
n=n+1
acof(n,1)=0.0
acof(n,2)=1.0
acof(n,3)=0.0
arhs(n)=v(101)
CALL LEQT1B(ACOF,N,1,1,IACOF,ARHS,1,IARHS,0,XL,IER)
C IACOF IS ACTUAL 1-ST DIMENSION OF THE ACOF ARRAY
C IARHS IS ACTUAL 1-ST DIMENSION OF THE ARHS ARRAY
C XL IS WORKING ARRAY WITH SIZE N*(NLC+1)
C WHERE N IS NUMBER OF EQUATIONS NLC NUMBER OF LOWER
C CODIAGONALS
C IF IER=129( .NE.0) MATRIX ACOF IS ALGORITHMICALLY SINGULAR
IF(IER.NE.0) CALL ERROUT
n=1
u(1)=arhs(n)
n=n+1
u(2)=arhs(n)
DO 25011 k=3,99,1
n=n+1
u(k)=arhs(n)
25011 CONTINUE
n=n+1
u(100)=arhs(n)
n=n+1
u(101)=arhs(n)
amer=0.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.0
arhs(n)=v(1)
n=n+1
alcd(n)=1.0
ad(n)=-2.0
aucd(n)=1.0
arhs(n)=v(3)-(2.0*v(2))+v(1)
DO 25014 k=3,99,1
n=n+1
alcd(n)=1.0
ad(n)=-2.0
aucd(n)=1.0
arhs(n)=v(k-1)+v(k+1)-(2.0*v(k))
25014 CONTINUE
n=n+1
alcd(n)=1.0
ad(n)=-2.0
aucd(n)=1.0
arhs(n)=v(101)-(2.0*v(100))+v(99)
n=n+1
alcd(n)=0.0
ad(n)=1.0
arhs(n)=v(101)
CALL DGTF(N,ALCD,AD,AUCD,AF,IPVT)
C N IS NUMBER OF EQUATIONS
C ALCD,AD,AUCD,AF,ARHS ARE ARRAYS OF DIMENSION AT LEAST (N)
C THESE ARRAYS ARE DOUBLE PRECISION TYPE
C FOR SINGLE PRECISION CHANGE DGTF TO SGTF AND DGTS TO SGTS
C IPVT IS INTEGER ARRAY OF DIMENSION AT LEAST (N+3)/8
CALL DGTS(N,ALCD,AD,AUCD,AF,IPVT,ARHS)
n=1
u(1)=arhs(n)
n=n+1
u(2)=arhs(n)
DO 25015 k=3,99,1
n=n+1
u(k)=arhs(n)
25015 CONTINUE
n=n+1
u(100)=arhs(n)
n=n+1
u(101)=arhs(n)
amer=0.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 17/05/91
*****************************
Partial Differential Equations
==============================
diff(v,t) - diff(w,x)*c = 0
- diff(v,x)*c + diff(w,t) = 0
0 interpolations are needed in t coordinate
Equation for v variable is integrated in half grid point
Equation for w variable is integrated in half grid point
0 interpolations are needed in x coordinate
Equation for v variable is integrated in 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
((w(n + 1,---------) - w(n + 1,---------) + w(n,---------)
2 2 2
2*j + 1
- w(n,---------))*c*ht + 2*v(n + 1,j)*hx - 2*v(n,j)*hx)/(2*ht*hx)
2
= 0
( - (v(n + 1,j + 1) - v(n + 1,j) + v(n,j + 1) - v(n,j))*c*ht
2*j + 1 2*j + 1
+ 2*w(n + 1,---------)*hx - 2*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
+ v(n,j)*c*ht + 2*w(n + 1,---------)*hx
2
2*j + 1
- 2*w(n,---------)*hx)/(2*ht*hx)
2
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*sin(ax) *c *ht *hx + hx )/(sin(ax) *c *ht
2 2 2 2 4
+ 2*sin(ax) *c *ht *hx + 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*sin(ax) *c *ht *hx + hx )/(sin(ax) *c *ht
2 2 2 2 4
+ 2*sin(ax) *c *ht *hx + hx )
pol:=hurw num pol;
2 2 2 2 2 2 2 2
pol := 4*lam *sin(ax) *c *ht *(sin(ax) *c *ht + hx )
2 2 2 2 2
+ 4*hx *(sin(ax) *c *ht + hx )
hurwitzp pol;
Necessary and sufficient conditions are:
2
hx
----------------- > 0
2 2 2
sin(ax) *c *ht
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
mat((-------------------------,---------------------),
2 3
hx hx
3 3 3
- 8*sin(ax) *c *ht *i
(------------------------,
3
hx
4 4 4 2 2 2 2 4
16*sin(ax) *c *ht - 4*sin(ax) *c *ht *hx + hx
--------------------------------------------------))
4
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
b := mat((-------------------------------------------,
hx
- i*sin(ax)*c*ht - i*sin(ay)*c*ht
-------------------,-------------------),
hx 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) + lam
2 2 2 2
*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
2 2 2
+ 4*cos(ay)*c*ht*hx + 7*c *ht - 8*c*ht*hx + 3*hx )
3 3 2 2
+ 4*cos(ax)*cos(ay)*c *ht - 3*cos(ax)*cos(ay)*c *ht *hx
3 3 2 2
- 4*cos(ax)*c *ht + 5*cos(ax)*c *ht *hx
2 3 3
- 2*cos(ax)*c*ht*hx - 4*cos(ay)*c *ht
2 2 2 3 3
+ 5*cos(ay)*c *ht *hx - 2*cos(ay)*c*ht*hx + 4*c *ht
2 2 2 3 3
- 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 2 2 2
pol := lam + lam *(4*s1 *cap + 4*s2 *cap - 3) + lam*(12*s1 *s2 *cap
2 2 2 2 2 2
+ 4*s1 *cap - 8*s1 *cap + 4*s2 *cap - 8*s2 *cap + 3)
2 2 3 2 2 2 2 2 2
+ 16*s1 *s2 *cap - 12*s1 *s2 *cap - 4*s1 *cap + 4*s1 *cap
2 2 2
- 4*s2 *cap + 4*s2 *cap - 1
clear cos ax,cos ay,sin ax,sin ay,cos ax2**2,cos ay2**2,sin ax2,sin ay2,
hx,hy;
pol:=hurw num pol;
3 2 2 3
pol := 16*lam *s1 *s2 *cap
2 2 2 2 2 2 2 2
+ 8*lam *cap *( - 6*s1 *s2 *cap + 3*s1 *s2 + s1 + s2 ) + 16
2 2 2 2 2 2 2
*lam*cap*(3*s1 *s2 *cap - 3*s1 *s2 *cap - s1 *cap + s1
2 2 2 2 3
- s2 *cap + s2 ) + 8*( - 2*s1 *s2 *cap
2 2 2 2 2 2 2 2
+ 3*s1 *s2 *cap + s1 *cap - 2*s1 *cap + s2 *cap
2
- 2*s2 *cap + 1)
hurwitzp pol;
If we denote:
2 2 3
(1) 16*s1 *s2 *cap > 0
2 2 2 2 2 2 2
(2) 8*cap *( - 6*s1 *s2 *cap + 3*s1 *s2 + s1 + s2 ) > 0
(3) 16*cap
2 2 2 2 2 2 2 2 2
*(3*s1 *s2 *cap - 3*s1 *s2 *cap - s1 *cap + s1 - s2 *cap + s2 )
> 0
2 2 3 2 2 2 2 2 2
(4) 8*( - 2*s1 *s2 *cap + 3*s1 *s2 *cap + s1 *cap - 2*s1 *cap
2 2 2
+ s2 *cap - 2*s2 *cap + 1) > 0
2 2 2 2 2 2 2
(5) 8*cap *( - 6*s1 *s2 *cap + 3*s1 *s2 + s1 + s2 ) > 0
3 4 4 3 4 4 2 4 4
(6) 128*cap *( - 16*s1 *s2 *cap + 24*s1 *s2 *cap - 9*s1 *s2 *cap
4 2 2 4 2 4 2
+ 8*s1 *s2 *cap - 10*s1 *s2 *cap + 3*s1 *s2
4 4 2 4 2 2 4
- s1 *cap + s1 + 8*s1 *s2 *cap - 10*s1 *s2 *cap
2 4 2 2 2 2 4
+ 3*s1 *s2 - 2*s1 *s2 *cap + s1 *s2 - s2 *cap
4
+ s2 ) > 0
3 6 6 6 6 6 5 6 6 4
(7) 1024*cap *(32*s1 *s2 *cap - 96*s1 *s2 *cap + 90*s1 *s2 *cap
6 6 3 6 4 5
- 27*s1 *s2 *cap - 32*s1 *s2 *cap
6 4 4 6 4 3
+ 100*s1 *s2 *cap - 93*s1 *s2 *cap
6 4 2 6 2 4
+ 27*s1 *s2 *cap + 10*s1 *s2 *cap
6 2 3 6 2 2 6 2
- 31*s1 *s2 *cap + 26*s1 *s2 *cap - 6*s1 *s2 *cap
6 3 6 2 6
- s1 *cap + 3*s1 *cap - 2*s1 *cap
4 6 5 4 6 4
- 32*s1 *s2 *cap + 100*s1 *s2 *cap
4 6 3 4 6 2
- 93*s1 *s2 *cap + 27*s1 *s2 *cap
4 4 4 4 4 3
+ 20*s1 *s2 *cap - 76*s1 *s2 *cap
4 4 2 4 4 4 2 3
+ 73*s1 *s2 *cap - 21*s1 *s2 *cap - 3*s1 *s2 *cap
4 2 2 4 2 4 2
+ 16*s1 *s2 *cap - 14*s1 *s2 *cap + 3*s1 *s2
4 4 2 6 4 2 6 3
- s1 *cap + s1 + 10*s1 *s2 *cap - 31*s1 *s2 *cap
2 6 2 2 6 2 4 3
+ 26*s1 *s2 *cap - 6*s1 *s2 *cap - 3*s1 *s2 *cap
2 4 2 2 4 2 4
+ 16*s1 *s2 *cap - 14*s1 *s2 *cap + 3*s1 *s2
2 2 2 2 6 3 6 2
- 2*s1 *s2 *cap + s1 *s2 - s2 *cap + 3*s2 *cap
6 4 4
- 2*s2 *cap - s2 *cap + s2 ) > 0
(c1) (1) AND (2) AND (4)
(c2) (1) AND (3) AND (4)
(d1) (5) AND (7)
(d2) (6)
Necessary and sufficient conditions are:
( (C1) OR (C2) ) AND ( (D1) OR (D2) )
cond
bt:=tcon b;
cos(ax)*c*ht + cos(ay)*c*ht - 2*c*ht + hx
bt := mat((-------------------------------------------,
hx
sin(ax)*c*ht*i sin(ay)*c*ht*i
----------------,----------------),
hx 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
- 4*cos(ay)*c *ht + 2*cos(ay)*c*ht*hx + 6*c *ht - 4*c*ht*hx
2 2
2 2 sin(ax)*c *ht *i*( - cos(ay) + 1)
+ hx )/hx ,-----------------------------------,
2
hx
2 2
sin(ay)*c *ht *i*( - cos(ax) + 1)
-----------------------------------),
2
hx
2 2
sin(ax)*c *ht *i*(cos(ay) - 1) 2 2
(--------------------------------,( - 2*cos(ax)*c *ht
2
hx
2 2 2 2
+ 2*cos(ax)*c*ht*hx + 2*c *ht - 2*c*ht*hx + hx )/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*cos(ay)*c *ht + 2*cos(ay)*c*ht*hx + 2*c *ht
2 2
- 2*c*ht*hx + hx )/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;
(fide 45915 2500)
End of Lisp run after 45.93+3.64 seconds