1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
|
/* ode1_lagrange.mac
Solution of Lagrange equation y = x f(y') + g(y')
References:
Daniel Zwillinger, Handbook of Differential Equations, 3rd ed
Academic Press, (1997), pp 328-331
Copyright (C) 2004 David Billinghurst
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
declare(method,special);
put('ode1_lagrange,001,'version)$
ode1_lagrange(eq,y,x):= block(
[eqn,ans1,ans:[],e,t,f,g,t2,s,p,q,%t],
ode_disp("In ode1_lagrange"),
/* Express equation as y = x*f(p)+g(p), where p=y' */
eqn:subst(p,'diff(y,x),eq),
if freeof(y,eqn) then return(false),
eqn:solve(eqn,y),
if eqn=[] then return(false),
/* There may be multiple solutions */
for e in eqn do (
ans1:block(
if lhs(e)#y then return(false),
if not(freeof(y,rhs(e))) then return(false),
t:rhs(e),
f:ratcoeff(t,x),
g:t-f*x,
if not(freeof(x,y,f)) then return(false),
if not(freeof(x,y,g)) then return(false),
if f=p then return(false),
/* Transform ode */
t2:'diff(x,p)=x*diff(f,p)/(p-f)+diff(g,p)/(p-f),
ode_disp2("ode1_lagrange: Transform successful and new equation is ",t2),
s:ode2(t2,x,p),
if (s=false) then return(false),
ode_disp("Transformed equation solved"),
method:'lagrange,
t:y=x*f+g,
/* Simple test to prevent failures */
ode_disp("Need to eliminate p from"),
ode_disp2(" Equation t: ",t),
ode_disp2(" Solution s: ",s),
/* Eliminate is only for polynomials.
Punt that it returns 1 when called inappropriately */
q:first(eliminate([t,s],[p])),
ode_disp(" and after eliminating p the result is"),
ode_disp2(" q: ",q),
if ( freeof(p,q) and (q#1) ) then (
q=0
)
else (
/* Return the parametric equation in %t=p */
[subst(%t,p,s),subst(%t,p,t)]
)
),
if ans1#false then ans:endcons(ans1,ans)
),
if ans=[] then
false
else
ans
)$
|
