1depend y,x;
2generic_function f(x,y);
3df(f(),x);
4df(f(x,y),x);
5df(f(x,x**3),x);
6df(f(x,z**3),x);
7df(a*f(x,y),x);
8dfp(a*f(x,y),x);
9df(f(x,y),x,2);
10df(dfp(f(x,y),x),x);
11df(dfp(f(x,x**3),x),x);
12
13% using a generic fucntion with commutative derivatives
14generic_function u(x,y);
15dfp_commute u(x,y);
16df(u(x,y),x,x);
17
18% explicitly declare 1st and second derivative commutative
19generic_function v(x,y);
20let dfp(v(~a,~b),{y,x}) => dfp(v(a,b),{x,y});
21df(v(),x,2);
22
23% substitute expressions for the arguments
24w:=df(f(),x,2);
25sub(x=0,y=x,w);
26
27% composite generic functions
28generic_function g(x,y);
29generic_function h(y,z);
30depend z,x;
31w:=df(g()*h(),x);
32sub(y=0,w);
33% substituting g*h for f in a partial derivative of f,
34% inheriting the arguments of f. Here no derivative of h
35% appears because h does not depend of x.
36sub(f=g*h,dfp(f(a,b),x));
37
38% indexes.
39
40% in the following total differential the partial
41% derivatives wrt i and j do not appear because i and
42% j do not depend of x.
43
44generic_function m(i,j,x,y);
45df(m(i,j,x,y),x);
46
47% computation with a differential equation.
48
49generic_function f(x,y);
50operator y;
51let df(y(~x),x) => f(x,y(x));
52
53% some derivatives
54
55df(y(x),x);
56df(y(x),x,2);
57df(y(x),x,3);
58sub(x=22,ws);
59
60% taylor expansion for y
61
62load_package taylor;
63taylor(y(x0+h),h,0,3);
64
65clear w;
66
67%------------------------ Runge Kutta -------------------------
68% computing Runge Kutta formulas for ODE systems Y'=F(x,y(x));
69% forms corresponding to Ralston Rabinowitz
70
71load_package taylor;
72operator alpha,betta,w,k;
73
74% s= order of Runge Kutta formula
75
76s:=3;
77
78generic_function f(x,y);
79operator y;
80
81% introduce ODE
82
83let df(y(~x),x)=>f(x,y(x));
84
85% formal series for solution
86
87y1_form := taylor(y(x0+h),h,0,s);
88
89% Runge-Kutta Ansatz:
90
91let alpha(1)=>0;
92
93for i:=1:s do
94   let k(i) => h*f(x0 + alpha(i)*h,
95                  y(x0) + for j:=1:(i-1) sum betta(i,j)*k(j));
96y1_ansatz:= y(x0) + for i:=1:s sum w(i)*k(i);
97
98y1_ansatz := taylor(y1_ansatz,h,0,s);
99
100% compute y1_form - y1_ans and collect coeffients of powers of h
101
102y1_diff := num(taylortostandard(y1_ansatz)-taylortostandard(y1_form))$
103cl := coeff(y1_diff,h);
104
105% f_forms: forms of f and its derivatives which occur in cl
106
107f_forms :=q := {f(x0,y(x0))}$
108for i:=1:(s-1) do
109  <<q:= for each r in q join {dfp(r,x),dfp(r,y)};
110    f_forms := append(f_forms,q);
111  >>;
112f_forms;
113
114% extract coefficients of the f_forms in cl
115
116sys := cl$
117for each fr in f_forms do
118  sys:=for each c in sys join coeff(c,fr);
119% and eliminate zeros
120sys := for each c in sys join if c neq 0 then {c} else {};
121
122end;
123