Artifact 36f5a003125227e23813741a4b74501fed7c67c55402bd67d7fddf7dd135b93d:
- File
r36/xlog/DFPART.LOG
— part of check-in
[152fb3bdbb]
at
2011-10-17 17:58:33
on branch master
— svn:eol-style, svn:executable and line endings for files
in historical/r36 treegit-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1480 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: schoepf@users.sourceforge.net, size: 12489) [annotate] [blame] [check-ins using] [more...]
REDUCE 3.6, 15-Jul-95, patched to 6 Mar 96 ... depend y,x; generic_function f(x,y); df(f(),x); df(y,x)*f (x,y) + f (x,y) y x df(f(x,y),x); df(y,x)*f (x,y) + f (x,y) y x df(f(x,x**3),x); 3 3 2 f (x,x ) + 3*f (x,x )*x x y df(f(x,z**3),x); 3 f (x,z ) x df(a*f(x,y),x); a*(df(y,x)*f (x,y) + f (x,y)) y x dfp(a*f(x,y),x); f (x,y)*a x df(f(x,y),x,2); 2 df(y,x,2)*f (x,y) + df(y,x) *f (x,y) + df(y,x)*f (x,y) + df(y,x)*f (x,y) y yy xy yx + f (x,y) xx df(dfp(f(x,y),x),x); df(y,x)*f (x,y) + f (x,y) xy xx df(dfp(f(x,x**3),x),x); 3 3 2 f (x,x ) + 3*f (x,x )*x xx xy % using a generic fucntion with commutative derivatives generic_function u(x,y); dfp_commute u(x,y); df(u(x,y),x,x); 2 df(y,x,2)*u (x,y) + df(y,x) *u (x,y) + 2*df(y,x)*u (x,y) + u (x,y) y yy xy xx % explicitly declare 1st and second derivative commutative generic_function v(x,y); let dfp(v(~a,~b),{y,x}) => dfp(v(a,b),{x,y}); df(v(),x,2); 2 df(y,x,2)*v (x,y) + df(y,x) *v (x,y) + 2*df(y,x)*v (x,y) + v (x,y) y yy xy xx % substitute expressions for the arguments w:=df(f(),x,2); 2 w := df(y,x,2)*f (x,y) + df(y,x) *f (x,y) + df(y,x)*f (x,y) + df(y,x)*f (x,y) y yy xy yx + f (x,y) xx sub(x=0,y=x,w); f (0,x) + f (0,x) + f (0,x) + f (0,x) xx xy yx yy % composite generic functions generic_function g(x,y); generic_function h(y,z); depend z,x; w:=df(g()*h(),x); w := df(y,x)*g (x,y)*h() + df(y,x)*h (y,z)*g() + df(z,x)*h (y,z)*g() + g (x,y)*h() y y z x sub(y=0,w); df(z,x)*h (0,z)*g(x,0) + g (x,0)*h(0,z) z x % substituting g*h for f in a partial derivative of f, % inheriting the arguments of f. Here no derivative of h % appears because h does not depend of x. sub(f=g*h,dfp(f(a,b),x)); g (a,b)*h(b,z) x % indexes. % in the following total differential the partial % derivatives wrt i and j do not appear because i and % j do not depend of x. generic_function m(i,j,x,y); df(m(i,j,x,y),x); df(y,x)*m (i,j,x,y) + m (i,j,x,y) y x % computation with a differential equation. generic_function f(x,y); operator y; let df(y(~x),x) => f(x,y(x)); % some derivatives df(y(x),x); f(x,y(x)) df(y(x),x,2); f (x,y(x)) + f (x,y(x))*f(x,y(x)) x y df(y(x),x,3); f (x,y(x)) + f (x,y(x))*f(x,y(x)) + f (x,y(x))*f (x,y(x)) xx xy x y 2 2 + f (x,y(x))*f(x,y(x)) + f (x,y(x))*f(x,y(x)) + f (x,y(x)) *f(x,y(x)) yx yy y sub(x=22,ws); f (22,y(22)) + f (22,y(22))*f(22,y(22)) + f (22,y(22))*f (22,y(22)) xx xy x y 2 + f (22,y(22))*f(22,y(22)) + f (22,y(22))*f(22,y(22)) yx yy 2 + f (22,y(22)) *f(22,y(22)) y % taylor expansion for y load_package taylor; taylor(y(x0+h),h,0,3); f (x0,y(x0)) + f (x0,y(x0))*f(x0,y(x0)) x y 2 y(x0) + f(x0,y(x0))*h + -----------------------------------------*h + ( 2 f (x0,y(x0)) + f (x0,y(x0))*f(x0,y(x0)) + f (x0,y(x0))*f (x0,y(x0)) xx xy x y 2 + f (x0,y(x0))*f(x0,y(x0)) + f (x0,y(x0))*f(x0,y(x0)) yx yy 2 3 4 + f (x0,y(x0)) *f(x0,y(x0)))/6*h + O(h ) y clear w; %------------------------ Runge Kutta ------------------------- % computing Runge Kutta formulas for ODE systems Y'=F(x,y(x)); % forms corresponding to Ralston Rabinowitz load_package taylor; operator alpha,beta,w,k; % s= order of Runge Kutta formula s:=3; s := 3 generic_function f(x,y); operator y; *** y already defined as operator % introduce ODE let df(y(~x),x)=>f(x,y(x)); % formal series for solution y1_form := taylor(y(x0+h),h,0,s); f (x0,y(x0)) + f (x0,y(x0))*f(x0,y(x0)) x y 2 y1_form := y(x0) + f(x0,y(x0))*h + -----------------------------------------*h 2 + (f (x0,y(x0)) + f (x0,y(x0))*f(x0,y(x0)) xx xy + f (x0,y(x0))*f (x0,y(x0)) + f (x0,y(x0))*f(x0,y(x0)) x y yx 2 2 3 + f (x0,y(x0))*f(x0,y(x0)) + f (x0,y(x0)) *f(x0,y(x0)))/6*h yy y 4 + O(h ) % Runge-Kutta Ansatz: let alpha(1)=>0; for i:=1:s do let k(i) => h*f(x0 + alpha(i)*h, y(x0) + for j:=1:(i-1) sum beta(i,j)*k(j)); y1_ansatz:= y(x0) + for i:=1:s sum w(i)*k(i); y1_ansatz := f(alpha(3)*h + x0, beta(3,2)*f(alpha(2)*h + x0,beta(2,1)*f(x0,y(x0))*h + y(x0))*h + beta(3,1)*f(x0,y(x0))*h + y(x0))*w(3)*h + f(alpha(2)*h + x0,beta(2,1)*f(x0,y(x0))*h + y(x0))*w(2)*h + f(x0,y(x0))*w(1)*h + y(x0) y1_ansatz := taylor(y1_ansatz,h,0,s); y1_ansatz := y(x0) + f(x0,y(x0))*(w(3) + w(2) + w(1))*h + ( alpha(3)*f (x0,y(x0))*w(3) + alpha(2)*f (x0,y(x0))*w(2) x x + beta(3,2)*f (x0,y(x0))*f(x0,y(x0))*w(3) y + beta(3,1)*f (x0,y(x0))*f(x0,y(x0))*w(3) y 2 + beta(2,1)*f (x0,y(x0))*f(x0,y(x0))*w(2))*h + ( y 2 alpha(3) *f (x0,y(x0))*w(3) xx + alpha(3)*beta(3,2)*f (x0,y(x0))*f(x0,y(x0))*w(3) xy + alpha(3)*beta(3,2)*f (x0,y(x0))*f(x0,y(x0))*w(3) yx + alpha(3)*beta(3,1)*f (x0,y(x0))*f(x0,y(x0))*w(3) xy + alpha(3)*beta(3,1)*f (x0,y(x0))*f(x0,y(x0))*w(3) yx 2 + alpha(2) *f (x0,y(x0))*w(2) xx + 2*alpha(2)*beta(3,2)*f (x0,y(x0))*f (x0,y(x0))*w(3) x y + alpha(2)*beta(2,1)*f (x0,y(x0))*f(x0,y(x0))*w(2) xy + alpha(2)*beta(2,1)*f (x0,y(x0))*f(x0,y(x0))*w(2) yx 2 2 + beta(3,2) *f (x0,y(x0))*f(x0,y(x0)) *w(3) yy 2 + 2*beta(3,2)*beta(3,1)*f (x0,y(x0))*f(x0,y(x0)) *w(3) yy 2 + 2*beta(3,2)*beta(2,1)*f (x0,y(x0)) *f(x0,y(x0))*w(3) y 2 2 + beta(3,1) *f (x0,y(x0))*f(x0,y(x0)) *w(3) yy 2 2 3 4 + beta(2,1) *f (x0,y(x0))*f(x0,y(x0)) *w(2))/2*h + O(h ) yy % compute y1_form - y1_ans and collect coeffients of powers of h y1_diff := num(taylortostandard(y1_ansatz)-taylortostandard(y1_form))$ cl := coeff(y1_diff,h); cl := {0, 6*f(x0,y(x0))*(w(3) + w(2) + w(1) - 1), 3*(2*alpha(3)*f (x0,y(x0))*w(3) + 2*alpha(2)*f (x0,y(x0))*w(2) x x + 2*beta(3,2)*f (x0,y(x0))*f(x0,y(x0))*w(3) y + 2*beta(3,1)*f (x0,y(x0))*f(x0,y(x0))*w(3) y + 2*beta(2,1)*f (x0,y(x0))*f(x0,y(x0))*w(2) - f (x0,y(x0)) y x - f (x0,y(x0))*f(x0,y(x0))), y 2 3*alpha(3) *f (x0,y(x0))*w(3) xx + 3*alpha(3)*beta(3,2)*f (x0,y(x0))*f(x0,y(x0))*w(3) xy + 3*alpha(3)*beta(3,2)*f (x0,y(x0))*f(x0,y(x0))*w(3) yx + 3*alpha(3)*beta(3,1)*f (x0,y(x0))*f(x0,y(x0))*w(3) xy + 3*alpha(3)*beta(3,1)*f (x0,y(x0))*f(x0,y(x0))*w(3) yx 2 + 3*alpha(2) *f (x0,y(x0))*w(2) xx + 6*alpha(2)*beta(3,2)*f (x0,y(x0))*f (x0,y(x0))*w(3) x y + 3*alpha(2)*beta(2,1)*f (x0,y(x0))*f(x0,y(x0))*w(2) xy + 3*alpha(2)*beta(2,1)*f (x0,y(x0))*f(x0,y(x0))*w(2) yx 2 2 + 3*beta(3,2) *f (x0,y(x0))*f(x0,y(x0)) *w(3) yy 2 + 6*beta(3,2)*beta(3,1)*f (x0,y(x0))*f(x0,y(x0)) *w(3) yy 2 + 6*beta(3,2)*beta(2,1)*f (x0,y(x0)) *f(x0,y(x0))*w(3) y 2 2 + 3*beta(3,1) *f (x0,y(x0))*f(x0,y(x0)) *w(3) yy 2 2 + 3*beta(2,1) *f (x0,y(x0))*f(x0,y(x0)) *w(2) - f (x0,y(x0)) yy xx - f (x0,y(x0))*f(x0,y(x0)) - f (x0,y(x0))*f (x0,y(x0)) xy x y 2 - f (x0,y(x0))*f(x0,y(x0)) - f (x0,y(x0))*f(x0,y(x0)) yx yy 2 - f (x0,y(x0)) *f(x0,y(x0))} y % f_forms: forms of f and its derivatives which occur in cl f_forms :=q := {f(x0,y(x0))}$ for i:=1:(s-1) do <<q:= for each r in q join {dfp(r,x),dfp(r,y)}; f_forms := append(f_forms,q); >>; f_forms; {f(x0,y(x0)), f (x0,y(x0)), x f (x0,y(x0)), y f (x0,y(x0)), xx f (x0,y(x0)), xy f (x0,y(x0)), yx f (x0,y(x0))} yy % extract coefficients of the f_forms in cl sys := cl$ for each fr in f_forms do sys:=for each c in sys join coeff(c,fr); % and eliminate zeros sys := for each c in sys join if c neq 0 then {c} else {}; sys := {6*(w(3) + w(2) + w(1) - 1), 3*(2*alpha(3)*w(3) + 2*alpha(2)*w(2) - 1), 3*(2*beta(3,2)*w(3) + 2*beta(3,1)*w(3) + 2*beta(2,1)*w(2) - 1), 2 2 3*alpha(3) *w(3) + 3*alpha(2) *w(2) - 1, 6*alpha(2)*beta(3,2)*w(3) - 1, 3*alpha(3)*beta(3,2)*w(3) + 3*alpha(3)*beta(3,1)*w(3) + 3*alpha(2)*beta(2,1)*w(2) - 1, 3*alpha(3)*beta(3,2)*w(3) + 3*alpha(3)*beta(3,1)*w(3) + 3*alpha(2)*beta(2,1)*w(2) - 1, 6*beta(3,2)*beta(2,1)*w(3) - 1, 2 2 3*beta(3,2) *w(3) + 6*beta(3,2)*beta(3,1)*w(3) + 3*beta(3,1) *w(3) 2 + 3*beta(2,1) *w(2) - 1} end; (TIME: dfpart 2390 2580)