Artifact 8143dd0449a63ec81594c10cc26a76fc22742da60decb856e51f109d8f653db6:
- Executable file
r36/XMPL/FIDE.TST
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 20390) [annotate] [blame] [check-ins using] [more...]
%*********************************************************************** %***** ***** %***** 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; a2:=div u; a3:=curl v; a4:=lapl v; a3:=2*a3+a4; a5:=lapl f; a1:=a1+div w; a1:=u.dyad((a,0,1),(1,b,3),(0,c,d)); a2:=vect(a,b,c); a1.a2; % Scalar product u.v; % Vector product u?v; % Dyadic u&v; % Directional derivative dirdf(u,v); 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; 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; 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; 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; 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; 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) ); % 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) ); % 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) ); % 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) ); % 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); % 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); % 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); %*********************************************************************** %***** ***** %***** 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; bb:=denotemat bb; factor lam; pol:=charpol bb; prdenot; 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; bb:=denotemat bb; pol:=charpol bb; prdenot; 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; bb:=denotemat bb; pol:=charpol bb; prdenot; 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; bb:=denotemat bb; pol:=charpol bb; prdenot; 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; pol:=charpol bb; 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; clear cos ax,cos ay,sin ax,sin ay,cos ax2**2,cos ay2**2,sin ax2,sin ay2, hx,hy; pol:=complexpol pol; pol1:=hurw pol; denotid cp; pol:=denotepol pol; prdenot; 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; 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; clear sin ax,cos ax,sin ay,cos ay,w,v,c,ht; pol:=charpol bb; denotid cp; pol:=denotepol pol; pol:=complexpol pol; denotid rp; pol:=denotepol pol; prdenot; cleardenot; clear aa,bb,pol; %*********************************************************************** %***** ***** %***** T e s t Examples --- Module H U R W P ***** %***** ***** %*********************************************************************** % Example H.1 x0:=lam-1; x1:=lam-(ar+i*ai); x2:=lam-(br+i*bi); x3:=lam-(cr+i*ci); hurwitzp x1; % Example H.2 x:=hurw(x0*x1); hurwitzp x; % Example H.3 x:=(x1*x2); hurwitzp x; % Example H.4 x:=(x1*x2*x3); hurwitzp x; 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; gentran <<literal tab!*,"dimension u(200),v(200),acof(200,3),arhs(200),xl(200,3)",cr!*$ dx:=0.05$ x:=0.1$ for i:=1:101 do <<v(i):=x**2/2$ x:=x+dx >> >>$ off period; gentran <<iacof:=200$ iarhs:=200 >>$ on period; genlinbandsol(1,1,{{u(1),u(1)=v(1)},{do,{k,2,100,1 },{u(k),u(k+1)- 2*u(k)+u(k-1)=v(k+1)-2*v(k)+v(k-1)}},{u(101),u(101)=v(101)}})$ gentran <<amer:=0.0$ arer:=0.0$ for i:=1:101 do <<am:=abs(u(i)-v(i))$ ar:=am/v(i)$ literal tab!*,"if(am.gt.amer) amer=am",cr!*$ literal tab!*,"if(ar.gt.arer) arer=ar",cr!* >>$ literal tab!*,"write(*,100)amer,arer",cr!*$ literal tab!*,"stop",cr!*$ literal "100 format(' max. abs. error = ',e12.2,", "' max. rel. error = ',e12.2)",cr!*$ literal tab!*,"end",cr!* >>$ on echo; %*********************************************************************** % Example L.2 on nag; off echo; gentran <<literal tab!*,"dimension u(200),v(200),acof(200,3),arhs(200),xl(200,3)",cr!*$ dx:=0.05$ x:=0.1$ for i:=1:101 do <<v(i):=x**2/2$ x:=x+dx >> >>$ off period; gentran <<iacof:=200$ iarhs:=200 >>$ on period; genlinbandsol(1,1,{{u(1),u(1)=v(1)},{do,{k,2,100,1 },{u(k),u(k+1)- 2*u(k)+u(k-1)=v(k+1)-2*v(k)+v(k-1)}},{u(101),u(101)=v(101)}})$ gentran <<amer:=0.0$ arer:=0.0$ for i:=1:101 do <<am:=abs(u(i)-v(i))$ ar:=am/v(i)$ literal tab!*,"if(am.gt.amer) amer=am",cr!*$ literal tab!*,"if(ar.gt.arer) arer=ar",cr!* >>$ literal tab!*,"write(*,100)amer,arer",cr!*$ literal tab!*,"stop",cr!*$ literal "100 format(' max. abs. error = ',e12.2,", "' max. rel. error = ',e12.2)",cr!*$ literal tab!*,"end",cr!* >>$ on echo; %*********************************************************************** % Example L.3 on imsl; off echo,nag; gentran <<literal tab!*,"dimension u(200),v(200),acof(200,3),arhs(200),xl(200,3)",cr!*$ dx:=0.05$ x:=0.1$ for i:=1:101 do <<v(i):=x**2/2$ x:=x+dx >> >>$ off period; gentran <<iacof:=200$ iarhs:=200 >>$ on period; genlinbandsol(1,1,{{u(1),u(1)=v(1)},{do,{k,2,100,1 },{u(k),u(k+1)- 2*u(k)+u(k-1)=v(k+1)-2*v(k)+v(k-1)}},{u(101),u(101)=v(101)}})$ gentran <<amer:=0.0$ arer:=0.0$ for i:=1:101 do <<am:=abs(u(i)-v(i))$ ar:=am/v(i)$ literal tab!*,"if(am.gt.amer) amer=am",cr!*$ literal tab!*,"if(ar.gt.arer) arer=ar",cr!* >>$ literal tab!*,"write(*,100)amer,arer",cr!*$ literal tab!*,"stop",cr!*$ literal "100 format(' max. abs. error = ',e12.2,", "' max. rel. error = ',e12.2)",cr!*$ literal tab!*,"end",cr!* >>$ on echo; %*********************************************************************** % Example L.4 on essl; off echo,imsl; gentran <<literal tab!*,"dimension u(200),v(200),acof(200,3),arhs(200),xl(200,3)",cr!*$ dx:=0.05$ x:=0.1$ for i:=1:101 do <<v(i):=x**2/2$ x:=x+dx >> >>$ off period; gentran <<iacof:=200$ iarhs:=200 >>$ on period; genlinbandsol(1,1,{{u(1),u(1)=v(1)},{do,{k,2,100,1 },{u(k),u(k+1)- 2*u(k)+u(k-1)=v(k+1)-2*v(k)+v(k-1)}},{u(101),u(101)=v(101)}})$ gentran <<amer:=0.0$ arer:=0.0$ for i:=1:101 do <<am:=abs(u(i)-v(i))$ ar:=am/v(i)$ literal tab!*,"if(am.gt.amer) amer=am",cr!*$ literal tab!*,"if(ar.gt.arer) arer=ar",cr!* >>$ literal tab!*,"write(*,100)amer,arer",cr!*$ literal tab!*,"stop",cr!*$ literal "100 format(' max. abs. error = ',e12.2,", "' max. rel. error = ',e12.2)",cr!*$ literal tab!*,"end",cr!* >>$ on echo; 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); on exp; center t=1/2; functions v,w; approx( aa(0,0)=aa(0,1)); center x=1/2; approx( aa(1,0)=aa(1,1)); 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); a(1,2):=aa(1,0); off prfourmat; b:=ampmat a; clear a,aa; factor lam; pol:=charpol b; pol:=troot1 pol; pol:=hurw num pol; hurwitzp pol; bt:=tcon b; bt*b; bt*b-b*bt; 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; clear a; factor lam; pol:=charpol b; pol:=hurw num pol; hurwitzp pol; bt:=tcon b; bt*b; bt*b-b*bt; 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); grid uniform,t,x,y; unfunc p,v1,v2; hy:=hx; off prfourmat; b:=ampmat a; pol:=charpol b; 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; 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; hurwitzp pol; bt:=tcon b; bt*b; bt*b-b*bt; clear a,b,bt,pol; %*********************************************************************** end;