#COPYRIGHT NOTICE: Copyright (c) 1997-1998 by Sasha Cyganowski. #Permission is #granted to anyone to use, modify, or redistribute this software freely, #subject to the following restrictions: 1. The author accepts no #responsibility for any consequences of this software and makes no guarantee #that this software is free of defects. 2. The origin of this software must #not be misrepresented, either by explicit claim or by omission. #3. This notice and the copyright must be included with all copies or #modified versions of this software. 4. This software may not be included or #redistributed as part of any package to be sold for profit unlesss the author #has given explicit written permission to do so. # Sasha Cyganowski, School of Computing and Mathematics, Deakin University, #Waurn Ponds, VIC, Australia. Internet: sash@deakin.edu.au # Telephone: +61 3 5227 1382. LJ:=proc(X,b,j) local k; sum('op(j,op(k,b))*diff(X,x[k])','k'=1..nops(b)); end: L0:=proc(X,a,b) local part1,part2,part3,k,j,l; part1:=diff(X,t); part2:=sum('a[k]*diff(X,x[k])','k'=1..nops(a)); part3:=1/2*sum('sum('sum('op(j,op(k,b))*op(j,op(l,b))*diff(X,x[k],x[l])','j'=1..nops(op(1,b)))','k'=1..nops(a))','l'=1..nops(a)); part1+part2+part3; end: WkTaylor2:=proc(a,b) local u,i,n,soln,Y,Delta,W,j1,j2,V; for i from 1 to nops(a) do soln[i]:=Y.i[n+1]=Y.i[n]+a[i]*Delta[n]+1/2*L0(a[i],a,b)*Delta[n]^2+sum('(op(j,op(i,b))+1/2*Delta[n]*(L0(op(j,op(i,b)),a,b)+LJ(a[i],b,j)))*Delta*W.j[n]','j'=1..nops(op(1,b)))+1/2*sum('sum('LJ(op(j2,op(i,b)),b,j1)*(Delta*W.j1[n]*Delta*W.j2[n]+V[(j1,j2)])','j1'=1..nops(op(1,b)))','j2'=1..nops(op(1,b))); for u from 1 to nops(a) do soln[i]:=subs(x[u]=Y.u[n],soln[i]); od; print(soln[i]); od; end: