#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. supp1:=proc(a,b,i) local soln,Y,Delta,n,j,u; soln[i]:=Y.i[n]+a[i]*Delta[n]+sum('op(j,op(i,b))*Delta*W.j[n]','j'=1..nops(op(1,b))); for u from 1 to nops(a) do soln[i]:=subs(x[u]=Y.u[n],soln[i]); od; end:\ supp2a:=proc(a,b,i,j) local soln,u,Y,n,Delta; soln[i]:=Y.i[n]+a[i]*Delta[n]+op(j,op(i,b))*sqrt(Delta[n]); for u from 1 to nops(a) do soln[i]:=subs(x[u]=Y.u[n],soln[i]); od; end: supp2b:=proc(a,b,i,j) local soln,u,Y,n,Delta; soln[i]:=Y.i[n]+a[i]*Delta[n]-op(j,op(i,b))*sqrt(Delta[n]); for u from 1 to nops(a) do soln[i]:=subs(x[u]=Y.u[n],soln[i]); od; end: supp3a:=proc(b,j,i) local soln,u,Y,n,Delta; soln[i]:=Y.i[n]+op(j,op(i,b))*sqrt(Delta[n]); for u from 1 to nops(b) do soln[i]:=subs(x[u]=Y.u[n],soln[i]); od; end: supp3b:=proc(b,j,i) local soln,u,Y,n,Delta; soln[i]:=Y.i[n]-op(j,op(i,b))*sqrt(Delta[n]); for u from 1 to nops(b) do soln[i]:=subs(x[u]=Y.u[n],soln[i]); od; end: supp4a:=proc(a,b,i,j) local tempb,k; tempb:=op(j,op(i,b)); for k from 1 to nops(a) do tempb:=subs(x[k]=supp2a(a,b,k,j),tempb); od; end: supp4b:=proc(a,b,i,j) local tempc,k; tempc:=op(j,op(i,b)); for k from 1 to nops(a) do tempc:=subs(x[k]=supp2b(a,b,k,j),tempc); od; end: supp5a:=proc(b,r,i,j) local tempd,k; tempd:=op(j,op(i,b)); for k from 1 to nops(b) do tempd:=subs(x[k]=supp3a(b,r,k),tempd); od; end: supp5b:=proc(b,r,i,j) local tempe,k; tempe:=op(j,op(i,b)); for k from 1 to nops(b) do tempe:=subs(x[k]=supp3b(b,r,k),tempe); od; end: expwk2:=proc(a,b) # Derivative free scheme local tempa,tempb,tempc,tempd,tempe,k,Y,Delta,W,j,i,n,V,soln,r,j1,r1,u,DeltaW; for i from 1 to nops(a) do tempa[i]:=op(i,a): for k from 1 to nops(a) do tempa[i]:=subs(x[k]=supp1(a,b,k),tempa[i]); od; soln[i]:=Y.i[n+1]=Y.i[n]+1/2*(tempa[i]+op(i,a))*Delta[n]+1/4*sum('((supp4a(a,b,i,j)+supp4b(a,b,i,j)+2*op(j,op(i,b)))*Delta*W.j[n]+(sum('(supp5a(b,r,i,j)+supp5b(b,r,i,j)-2*op(j,op(i,b)))*Delta*W.j[n]','r'=1..j-1)+sum('(supp5a(b,r,i,j)+supp5b(b,r,i,j)-2*op(j,op(i,b)))*Delta*W.j[n]','r'=j+1..nops(op(1,b)))))','j'=1..nops(op(1,b)))+1/4*sum('((supp4a(a,b,i,j1)-supp4b(a,b,i,j1))*((Delta*W.j1[n])^2-Delta[n])+(sum('(supp5a(b,r1,i,j1)-supp5b(b,r1,i,j1))*(DeltaW.j1[n]*(DeltaW.r1[n])+V[(j1,r1)])','r1'=1..j1-1)+sum('(supp5a(b,r1,i,j1)-supp5b(b,r1,i,j1))*(DeltaW.j1[n]*(DeltaW.r1[n])+V[(j1,r1)])','r1'=j1+1..nops(op(1,b)))))*(Delta[n])^(-1/2)','j1'=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: impwk2:=proc(a,b) # Derivative free scheme local tempa,k,Y,Delta,W,j,i,n,V,soln,r,j1,r1,u; for i from 1 to nops(a) do tempa[i]:=op(i,a): for k from 1 to nops(a) do tempa[i]:=subs(x[k]=supp1(a,b,k),tempa[i]); od; soln[i]:=Y.i[n+1]=Y.i[n]+1/2*(tempa[i]+op(i,a))*Delta[n]+1/4*sum('(supp4a(a,b,i,j)+supp4b(a,b,i,j)+2*op(j,op(i,b))+(sum('(supp5a(b,r,i,j)+supp5b(b,r,i,j)-2*op(j,op(i,b)))*(Delta[n])^(-1/2)','r'=1..j-1)+sum('(supp5a(b,r,i,j)+supp5b(b,r,i,j)-2*op(j,op(i,b)))*(Delta[n])^(-1/2)','r'=j+1..nops(op(1,b)))))*Delta*W.j[n]','j'=1..nops(op(1,b)))+1/4*sum('((supp4a(a,b,i,j1)-supp4b(a,b,i,j1))*((Delta*W.j1[n])^2-Delta[n])+(sum('(supp5a(b,r1,i,j1)-supp5b(b,r1,i,j1))*(DeltaW.j1[n]*(DeltaW.r1[n])+V[(j1,r1)])','r1'=1..j1-1)+sum('(supp5a(b,r1,i,j1)-supp5b(b,r1,i,j1))*(DeltaW.j1[n]*(DeltaW.r1[n])+V[(j1,r1)])','r1'=j1+1..nops(op(1,b)))))*(Delta[n])^(-1/2)','j1'=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: