6.5.2 Source code
#--------------------------------------------------------------------
#FILE NAME : KOV.mpl
#
# Copyright (C) 2022, Nasser M. Abbasi. All rights reserved
# email: nma@12000.org
#
# Free software to use and modify in anyway as long as the above
# copyright notice remains attached in the file
#
# Change history
#-----------------
# Oct 27, 2022. Initial version
#
# Note that the latest version and any updates can alaways be obtained
# from the author web site at
# http://12000.org/my_notes/my_paper_on_kovacic/paper.htm
#
# Any problems found in the software please report so I can correct.
# This implementation was done using Maple 2022.2 on windows 10.
#
#--------------------------------------------------------------------
#--------------------------------------------------------------------
# An Object oriented implementation of the Kovacic algortithm
# using Maple 2022.
#
# based on original Kovacic paper description of the algorithm.
#
# This file contains two modules. The first is called kovacic_class
# used to create kovacic object. The second is kovacic_tester module
# used for unit testing the kovacic_class module
#
# To use this file just do
#
# read "KOV.mpl"
# ode := diff(y(x),x$2)=2/x^2*y(x)
# o := Object(kovacic_class,ode,y(x));
# sol := o:-dsolve();
#
# Make sure to set the currentdir() in Maple correctly in order to find
# where you downloaded the file "KOV.mpl" to.
#--------------------------------------------------------------------
#--------------------------------------------------------------------
# This class is used to solve an ode using Kovacic algorithm.
#
# Please see documantion section in the paper for additional
# information on using this class.
#--------------------------------------------------------------------
kovacic_class :=module()
option object;
#class to hold one entry in the gamma set for case 1
local case_one_gamma_entry := module()
option object;
export pole_location := 0;
export pole_order := 0;
export sqrt_r := 0;
export alpha_plus := 0;
export alpha_minus := 0;
export b:=0;
end module;
#class to hold O_infinity information for case 1
local case_one_O_inf := module()
option object;
export sqrt_r_inf := 0;
export alpha_plus_inf := 0;
export alpha_minus_inf := 0;
export a := 0;
export b := 0;
end module;
#class to hold one entry in the gamma set for case 2,3
local case_2_and_3_gamma_entry:=module()
option object;
export pole_location := 0;
export pole_order := 0;
export Ec::set := {};
export b := 0;
end module;
#--------------------------------------------------------------------
# PRIVATE variables for the kovacic class
# only methods inside this class can access these
#--------------------------------------------------------------------
local original_ode;
#coefficients of original linear ode A y''+ B y' + C y =0
local A,B,C;
#original ode dependent and independent variables
local y::symbol,x::symbol;
local modified_ode; #this is the z''=r*z ode;
local z::symbol; #dependent variable for the r_ode z''(x)=r*z(x)
local r,s,t; #where r=s/t;
local O_inf; #order of r at infinity. degree(s)-degree(t)
# poles or r. format is [ [pole,order],[pole,order],...]
# if no poles, for example r=x, then list is empty []
# this means pole order zero is not in the list, since no pole.
local poles_list::list := [];
#contains all possible cases detected. Hence [1] or [1,2] or [1,2,3]
#or remains empty [] for case 4.
local list_of_possible_cases::list := [];
local case_used_to_solve::integer:=-1; #this will have case used: 1,2 or 3
#this will have n=4,6 or 12 used for case 3 only
local n_used_for_case_3::integer:=-1;
#--------------------------------------------------------------------
# CONSTRUCTOR
#
# the input is the ode itself as first argument, and the
# dependent variable, y(x) for example, as second argument.
#--------------------------------------------------------------------
export ModuleCopy::static:= proc( _self,
proto::kovacic_class,
ode::`=`,func::function(name) ,$)
local A,B,C;
local x,y,r,z::nothing,s,t;
x,y,A,B,C := parse_and_validate_ode(ode,func);
_self:-original_ode := ode;
#force r to be relatively prime ratio of 2 polynomials
r := normal( (2*diff(B,x)*A-2*B*diff(A,x)+B^2-4*A*C)/(4*A^2));
if not type(r,'ratpoly'(anything,x)) then
ERROR("r= ", r ," is not polynomial or rational function in ",x);
fi;
s := numer(r);
t := denom(r);
#save all findings into the object private data
#so it can be used later
_self:-modified_ode := diff(z(x),x$2) = r*z(x);
_self:-O_inf := degree(t,x)-degree(s,x);
_self:-r := r;
_self:-x := x;
_self:-y := y;
_self:-z := z;
_self:-s := s;
_self:-t := t;
_self:-C := C;
_self:-B := B;
_self:-A := A;
#determine all poles and order
_self:-generate_poles_and_order_list();
#determine all possible kovacic cases
_self:-find_possible_cases();
return NULL;
end proc;
#--------------------------------------------------------------------
# This module private function is called by constructor to parse and
# validate the ode
#--------------------------------------------------------------------
local parse_and_validate_ode:=proc(ode::`=`,func::function(name),$)
local x,y,A,B,C,L::list;
local item,dep_variables_found;
if nops(func)<>1 then
ERROR("dependent variable ",func," has more than one argument");
fi;
y := op(0,func);
x := op(1,func);
#basic verification
if not has(ode,y) then ERROR("Supplied ",ode," has no ",y); fi;
if not has(ode,x) then ERROR("Supplied ", ode," has no ",x); fi;
if not has(ode,func) then ERROR("Supplied ",ode," has no ",func); fi;
#check it is second order
if PDEtools:-difforder(ode)<>2 then
ERROR("Only second order ode's can be used in Kocacic algorithm. ");
fi;
#check all dependent variables in ode which is y, match given func y(x)
try
dep_variables_found := PDEtools:-Library:-GetDepVars([y],ode);
catch:
ERROR(lastexception);
end try;
#o over dep_variables_found and check the
#independent variable is same as x i.e. ode can be y'(z)+y(z)=0 but
#function is y(x).
for item in dep_variables_found do
if not type(item,function) then
ERROR("Parse error. Expected ",func," found ",item," in", ode);
else
if op(1,item) <> x then
ERROR("Parse error. Expected ",func," found ",item," in",ode);
fi;
if nops(item)<>1 then
ERROR("Parse error. One argument allowed in ",
func," found ", item," in " , ode);
fi;
fi;
od;
#now go over all indents in ode and check that y shows as y(x)
#and not as just y as the PDEtools:-Library:-GetDepVars([_self:-y],ode)
#code above does not detect this. i.e. it does not check y'(x)+y=0
if numelems(indets(ode,identical(y))) > 0 then
ERROR("Parsing error, Can not have ",y," with no argument inside ",ode);
fi;
#check ode is linear in y(x)
if not has(DEtools:-odeadvisor(ode,func,['linear']),'_linear') then
ERROR("Only linear ode's can be used in Kocacic algorithm. ");
fi;
#extract coefficients of the ode
L:=DEtools:-convertAlg(ode,func); #this only works on linear ode's
if L[2]<>0 then ERROR("Not homogeneous ode"); fi;
#Finished parsing the ode. A y''+B y' + C y = 0
C := L[1,1];
B := L[1,2];
A := L[1,3];
if not type(A,'ratpoly'(anything,x)) or not type(B,'ratpoly'(anything,x))
or not type(C,'ratpoly'(anything,x)) then
error "ode coefficients are not rational functions of ",x;
fi;
return x,y,A,B,C;
end proc;
#--------------------------------------------------------------------
# main API. Called to solve the ode using Kovacic algorithm.
# returns the solution in the form of y(x)=.... or if no solution
# is found returns FAIL.
# see user guide how to use.
#--------------------------------------------------------------------
export dsolve::static:=proc(_self,$)
local sol:=FAIL, current_case_number::posint;
if nops(_self:-list_of_possible_cases) = 0 then
return FAIL;
fi;
#keep trying all possible cases, starting from case 1 to case 3.
#until one case succeed or all are tried
for current_case_number in _self:-list_of_possible_cases do
sol := _self:-dsolve_case(current_case_number);
if sol<>FAIL then
_self:-case_used_to_solve:=current_case_number;
return sol;
fi;
od;
return FAIL;
end proc;
#--------------------------------------------------------------------
# Called to solve the ode using a specific case number
# made public to allow user to solve using specific case
#--------------------------------------------------------------------
export dsolve_case::static:=proc(_self,case_number::posint,$)
local sol;
if case_number>3 then ERROR("Only case number 1,2,3 are allowed"); fi;
if nops(_self:-list_of_possible_cases)=0 then
ERROR("No possible cases detected for this ode");
fi;
if not member(case_number,_self:-list_of_possible_cases) then
ERROR("Case ", case_number,
" not one of possible cases ", _self:-list_of_possible_cases);
fi;
if case_number = 1 then
sol:= _self:-solve_case_1();
elif case_number = 2 then
sol:= _self:-solve_case_2();
else
sol:= _self:-solve_case_3();
fi;
_self:-case_used_to_solve := case_number;
return sol;
end proc;
#--------------------------------------------------------------------
# returns back the z''(x) = r(x) z(x) ode, which is the one
# actually solved by Kovacic algorithm
#--------------------------------------------------------------------
export get_z_ode::static:=proc(_self,$)
local x := _self:-x;
local r := _self:-r;
local z := _self:-z;
if _self:-s = 0 then
return diff(z(x),x$2) = 0;
else
return diff(z(x),x$2) = numer(r)%/denom(r) * z(x);
fi;
end proc;
#--------------------------------------------------------------------
# returns back the r term in the z''(x) = r(x) z(x) ode
#--------------------------------------------------------------------
export get_r::static:=proc(_self,$)
if _self:-s = 0 then
return 0;
else
return numer(_self:-r)%/denom(_self:-r);
fi;
end proc;
#--------------------------------------------------------------------
# returns s, where r = s/t and z'' = r z(t)
#--------------------------------------------------------------------
export get_s::static:=proc(_self,$)
return _self:-s;
end proc;
#--------------------------------------------------------------------
# returns t, where r = s/t and z'' = r z(t)
#--------------------------------------------------------------------
export get_t::static:=proc(_self,$)
return _self:-t;
end proc;
#--------------------------------------------------------------------
# Returns back the original ode ( A y''+ B y' + C y = 0 )
#--------------------------------------------------------------------
export get_y_ode::static:=proc(_self,$)
return _self:-original_ode;
end proc;
#--------------------------------------------------------------------
# returns list of poles of r, where z''(x) = r z(x)
# The list has format [ [pole location,order], ...]
#--------------------------------------------------------------------
export get_poles::static:=proc(_self,$)
return _self:-poles_list;
end proc;
#--------------------------------------------------------------------
# returns O_infinity of r, where z''(x) = r z(x)
#--------------------------------------------------------------------
export get_order_at_infinity::static:=proc(_self,$)
return _self:-O_inf;
end proc;
#--------------------------------------------------------------------
# returns list of possible kovacic cases possible.
#--------------------------------------------------------------------
export get_possible_cases::static:=proc(_self,$)::list;
return _self:-list_of_possible_cases;
end proc;
#--------------------------------------------------------------------
# returns actual case number used when solving ode.
# can be 1,2 or 3
# if no cases applicable, then -1 is returned
#--------------------------------------------------------------------
export get_case_used::static:=proc(_self,$)::integer;
return _self:-case_used_to_solve;
end proc;
#--------------------------------------------------------------------
# returns n used for case 3. Either 4,6, or 12.
# if not case 3 used, then -1 is returned.
#--------------------------------------------------------------------
export get_n_case_3::static:=proc(_self,$)::integer;
return _self:-n_used_for_case_3;
end proc;
#--------------------------------------------------------------------
# All functions below are private
#--------------------------------------------------------------------
#--------------------------------------------------------------------
#This proc find all possible kovacic cases. These can be 1,2 or 3.
#if none of these found, then empty list is returned, which is
#case 4 in the paper. For example, if case 1 is only possible,
#then [1] is returned. if case 1 and 2 are possible, then [1,2]
#is returned, if all three cases possible, then [1,2,3] is returned.
#If no cases possible then [] returned.
#--------------------------------------------------------------------
local find_possible_cases::static:=proc(_self,$)
local L::list := [];
local poles_order := convert(_self:-poles_list[..,2],set);
local T::set;
#check for case 1
T := select(Z-> Z>2,poles_order);
if nops( select(Z->type(Z,odd),T) )=0 and
(_self:-O_inf<0 and type(_self:-O_inf,even)) or
_self:-O_inf=0 or _self:-O_inf>1 then
L:= [1];
fi;
if nops(poles_order)>0 then #must have at least one pole for 2,3 cases
if nops(poles_order)<>0 then #can not have pole order 0, i.e. no poles
T:=select(Z-> Z>2,poles_order);
# r must have at least one pole that is either odd order greater
#than 2 or else has order 2
if nops( select(Z->type(Z,odd),T) )<>0
or member(2,poles_order) then
L:= [ op(L),2 ];
fi;
#check for case 3
if not member(0,poles_order)
and nops(select(Z-> Z>2,poles_order))=0
and _self:-O_inf>=2 then
L:= [ op(L),3 ];
fi;
fi;
fi;
_self:-list_of_possible_cases := L;
return NULL;
end proc;
#--------------------------------------------------------------------
#called to find all poles of r and the order
#of each pole. Uses Maple's sqrfree.
#--------------------------------------------------------------------
local generate_poles_and_order_list::static:=proc(_self,$)
local L::list := [];
local sol::list;
local current_sol;
local poles_list::list;
local current_pole;
local r := _self:-r, x := _self:-x;
poles_list := sqrfree(denom(r),x);
poles_list := poles_list[2,..]; #we do not need the overall factor
if nops(poles_list) = 0 then
_self:-poles_list:=[];
else
for current_pole in poles_list do
sol := solve(current_pole[1]=0,[x]);
sol := ListTools:-Flatten(sol);
for current_sol in sol do
L := [op(L), [rhs(current_sol) ,current_pole[2] ] ];
od;
od;
_self:-poles_list := L;
fi;
return NULL;
end proc:
#--------------------------------------------------------------------
#
# C A S E O N E I M P L E M E N T A T I O N
#
# returns ode solution using case 1, or FAIL is no solution exist
#--------------------------------------------------------------------
local solve_case_1::static:=proc(_self,$)
local O_infinity_set::kovacic_class:-case_one_O_inf;
local gamma_set::set(kovacic_class:-case_one_gamma_entry):={};
gamma_set, O_infinity_set := _self:-case_1_step_1();
return _self:-case_1_step_2(gamma_set, O_infinity_set);
end proc;
#--------------------------------------------------------------------
# called from _self:-solve_case_1()
#--------------------------------------------------------------------
local case_1_step_1::static:=proc(_self,$)::
set(kovacic_class:-case_one_gamma_entry),
kovacic_class:-case_one_O_inf;
local current_pole;
local e::kovacic_class:-case_one_gamma_entry;
local o::kovacic_class:-case_one_O_inf;
local b,a,v,i::integer;
local b_coeff_in_r,b_coeff_in_r_inf_square;
local x := _self:-x, r := _self:-r;
local N::integer;
local laurent_c;
local current_term;
local b_from_r,b_from_laurent_series;
#this contains all information generated for each pole of r
#this is what is called the set GAMMA in the paper and
#in the diagram of algorithm above.
local gamma_set::set(kovacic_class:-case_one_gamma_entry) := {};
if nops(_self:-poles_list) = 0 then
gamma_set := {};
else
for current_pole in _self:-poles_list do
e := Object(kovacic_class:-case_one_gamma_entry);
e:-pole_location := current_pole[1];
e:-pole_order := current_pole[2];
if e:-pole_order =1 then
e:-sqrt_r := 0;
e:-alpha_plus := 1;
e:-alpha_minus := 1;
gamma_set := { op(gamma_set), e };
elif e:-pole_order = 2 then
e:-sqrt_r := 0;
e:-b := _self:-b_partial_fraction(r,x,e:-pole_location,2);
e:-alpha_plus := 1/2+1/2*sqrt(1+4*e:-b);
e:-alpha_minus := 1/2-1/2*sqrt(1+4*e:-b);;
gamma_set := { op(gamma_set), e };
else
v := e:-pole_order/2;
e:-sqrt_r := 0;
for N from 2 to v do
laurent_c := _self:-laurent_coeff(sqrt(r),x,
e:-pole_location,v,N);
current_term := laurent_c/(x-e:-pole_location)^N;
e:-sqrt_r := e:-sqrt_r + current_term;
if N = v then
a := laurent_c;
fi;
od;
b_from_r :=_self:-b_partial_fraction(r,x,e:-pole_location,v+1);
b_from_laurent_series :=_self:-laurent_coeff(sqrt(r),x,
e:-pole_location,v,v+1);
e:-b := b_from_r - b_from_laurent_series;
e:-alpha_plus := 1/2*((e:-b)/a + v);
e:-alpha_minus := 1/2*(-(e:-b)/a + v);
gamma_set := { op(gamma_set), e };
fi;
od;
fi;
o := Object(case_one_O_inf);
if _self:-O_inf > 2 then
o:-sqrt_r_inf := 0;
o:-alpha_plus_inf := 0;
o:-alpha_minus_inf := 1;
elif _self:-O_inf=2 then
o:-sqrt_r_inf :=0;
o:-b := lcoeff(_self:-s) / lcoeff(_self:-t);
b := radsimp((1+4*o:-b)^(1/2));
o:-alpha_plus_inf := 1/2 + 1/2*b;
o:-alpha_minus_inf := 1/2 - 1/2*b;
else #order at infinity -2*v<= 0 which must be even
v := (-_self:-O_inf) / 2;
o:-sqrt_r_inf :=0;
for i from 0 to v do
laurent_c := _self:-laurent_coeff(sqrt(r),x,infinity,v,i);
o:-sqrt_r_inf := o:-sqrt_r_inf + laurent_c*x^i;
if i = v then
o:-a := laurent_c;
fi;
od;
b_coeff_in_r_inf_square := _self:-laurent_coeff(
o:-sqrt_r_inf^2,x,0,v,v-1);
b_coeff_in_r := _self:-get_coefficient_of_r(v-1);
o:-b := b_coeff_in_r - b_coeff_in_r_inf_square;
o:-alpha_plus_inf := 1/2*( (o:-b)/(o:-a) - v);
o:-alpha_minus_inf := 1/2*( -(o:-b)/(o:-a) - v);
fi;
return gamma_set,o;
end proc;
#--------------------------------------------------------------------
# Finds b coefficient in r for the case one only when v<=0
# using long division
#--------------------------------------------------------------------
local get_coefficient_of_r::static:=proc(_self,the_degree,$)
local r := _self:-r;
local x := _self:-x;
local c,R,t;
try
c := coeff(r,x,the_degree);
catch:
R := rem(numer(r),denom(r),x);
t := denom(r);
c := lcoeff(R)/lcoeff(t);
end try;
return c;
end proc;
#--------------------------------------------------------------------
# called from _self:-solve_case_1()
# This determines the set of d non-negative integers and
# corresponding w for each d.
#--------------------------------------------------------------------
local case_1_step_2::static:=proc(_self,
gamma_set::set(kovacic_class:-case_one_gamma_entry),
O_infinity::kovacic_class:-case_one_O_inf,$)
local d,N,K,number_of_poles::integer,current_alpha_infinity;
local item::list;
local the_sign::integer,sign_list::list;
local w;
local r_solution,y_solution;
local x := _self:-x;
#this will contains all data found for any nonnegative d. Each
#entry will be a list of this form
# [d, w]
local good_d_found := Array(1..0);
local B::Matrix; #this will contain good_d_found but as matrix
number_of_poles := nops(gamma_set);
sign_list := combinat:-permute([1$number_of_poles,-1$number_of_poles],
number_of_poles);
for K,current_alpha_infinity in
[O_infinity:-alpha_plus_inf,O_infinity:-alpha_minus_inf] do
for item in sign_list do
d := 0;
if number_of_poles>0 then
for N,the_sign in item do
if the_sign = -1 then
d := d-gamma_set[N]:-alpha_minus;
else
d := d-gamma_set[N]:-alpha_plus;
fi;
od;
fi;
d := simplify(d + current_alpha_infinity);
if type(d,'integer') and d >= 0 then
w := 0;
for N,the_sign in item do
if number_of_poles>0 then
if the_sign = -1 then
w := w +(-1)*gamma_set[N]:-sqrt_r +
(gamma_set[N]:-alpha_minus)/
(x - gamma_set[N]:-pole_location);
else
w := w + gamma_set[N]:-sqrt_r +
( gamma_set[N]:-alpha_plus)/
(x - gamma_set[N]:-pole_location);
fi;
fi;
od;
#this to get the sign in according to paper. First + then -
if K = 1 then
w := w + O_infinity:-sqrt_r_inf;
else
w := w - O_infinity:-sqrt_r_inf;
fi;
good_d_found ,= [ d, w];
fi;
od;
od;
if numelems(good_d_found) = 0 then return FAIL; fi;
#now convert the array to Matrix, and sort on d, so we
#start will the smallest
#degree d, which is the first column, as that will be most efficient.
#convert to set first, to remove any possible duplicat entries
#then convert to Matrix
convert(good_d_found,set);
B := convert(convert(%,list),Matrix);
B := B[sort(B[.., 1], 'output'= 'permutation')];
for N from 1 to LinearAlgebra:-RowDimension(B) do
r_solution := _self:-case_1_step_3(B[N,1], B[N,2]); #(d,w)
if r_solution <> FAIL then
y_solution := _self:-build_y_solution_from_r_solution(r_solution);
return y_solution;
fi;
od;
return FAIL;
end proc;
#--------------------------------------------------------------------
# called from _self:-case_1_step_2()
#--------------------------------------------------------------------
local case_1_step_3::static:=proc(_self,d::nonnegative,w,$)
local x := _self:-x;
local r := _self:-r;
local p;
local i::integer;
local a::nothing;
local eq;
local coeff_sol;
local tmp,W;
local final_result;
p := x^d;
for i from d-1 by -1 to 0 do
p := p + a[i] * x^i;
od;
#using original kovacis method. Not Smith.
eq := simplify( diff(p,x$2)+2*w*diff(p,x)+(diff(w,x)+w^2-r)*p) = 0;
if d = 0 then
if not evalb(eq) then return FAIL; fi;
else
#solve for coefficients
try
coeff_sol := timelimit(30,solve(
identity(eq,x), [seq(a[i],i=0..d-1)]));
catch:
return FAIL;
end try;
if nops(coeff_sol) = 0 then return FAIL; fi;
tmp := map(evalb,coeff_sol[1]);
if has(tmp,true) then return FAIL; fi;
p := eval(p,coeff_sol[1]); #to force a[i] solutions to update
fi;
W := diff(p,x)/p + w;
W := radsimp(W);
tmp := diff(W,x)+W^2;
if evalb( tmp = r) or is(tmp = r) or simplify(tmp-r)=0 then #can be used
try
tmp := int(w, x);
catch:
return FAIL;
end try;
if has(tmp,int) then return FAIL; fi;
final_result := simplify(p*exp(tmp));
if has(final_result,signum) or has(final_result,csgn) then
final_result := p*exp(tmp);
fi;
return final_result;
else
return FAIL;
fi;
end proc;
#--------------------------------------------------------------------
#
# C A S E T W O I M P L E M E N T A T I O N
#
# returns ode solution using case 2, or FAIL is no solution exist
#--------------------------------------------------------------------
local solve_case_2::static:=proc(_self,$)
local E_inf::set;
local gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry):={};
gamma_set, E_inf := _self:-case_2_step_1();
return _self:-case_2_step_2(gamma_set, E_inf );
end proc;
#--------------------------------------------------------------------
#
#--------------------------------------------------------------------
local case_2_step_1::static:=proc(_self,$)::
set(kovacic_class:-case_2_and_3_gamma_entry),set;
local current_pole;
local e::kovacic_class:-case_2_and_3_gamma_entry;
local E_inf::set;
local b;
local x := _self:-x;
local r := _self:-r;
#this contains all information generated for each pole of r
#this is what is called the set GAMMA in the paper and in the diagram of
#algorithm above.
local gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry) := {};
for current_pole in _self:-poles_list do
e := Object(kovacic_class:-case_2_and_3_gamma_entry);
e:-pole_location := current_pole[1];
e:-pole_order := current_pole[2];
if e:-pole_order = 1 then
e:-Ec := {4};
gamma_set := { op(gamma_set), e };
elif e:-pole_order = 2 then
e:-b := _self:-b_partial_fraction(r,x,e:-pole_location,2);
e:-Ec := {2,2+2*sqrt(1+4* e:-b),2-2*sqrt(1+4* e:-b)};
e:-Ec := select(z->type(z,integer),e:-Ec);
gamma_set := { op(gamma_set), e };
else
e:-Ec := {e:-pole_order};
gamma_set := { op(gamma_set), e };
fi;
od;
if _self:-O_inf>2 then
E_inf := {0,2,4};
elif _self:-O_inf=2 then
b := lcoeff(_self:-s) / lcoeff(_self:-t);
E_inf := {2,2+2*sqrt(1+4*b),2-2*sqrt(1+4*b)};
E_inf := select(z->type(z,integer),E_inf);
else #order at infinity v< 2
E_inf := {_self:-O_inf};
fi;
return gamma_set,E_inf;
end proc;
#--------------------------------------------------------------------
# called from _self:-solve_case_2()
# This determines the set of d non-negative integers and
# corresponding w for each d.
#--------------------------------------------------------------------
local case_2_step_2::static:=proc(_self,
gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry),
E_inf::set,
$)
local L::list := [];
local item,current_E_inf;
local d,N,ee,theta;
local x := _self:-x;
local r_solution,y_solution;
#this will contains all data found for any nonnegative d. Each
#entry will be a list of this form
# [d, theta]
#so if we obtain say 3 values of d that are nonnegative,
#there will be 3 such lists in this array
local good_d_found := Array(1..0);
local B::Matrix; #this will contain good_d_found as matrix
for item in gamma_set do
L := [ op(L), convert(item:-Ec,list) ];
od;
#now find all possible tuples
if nops(L)>1 then
L := kovacic_class:-cartProdSeq(op(L));
else
L := L[1];
fi;
for current_E_inf in E_inf do
for item in L do
d := 1/2*( current_E_inf - add(item));
if type(d,'integer') and d >= 0 then
theta :=0;
for N,ee in item do
theta := theta + ee/(x-gamma_set[N]:-pole_location);
od;
theta := 1/2*theta;
good_d_found ,= [ d, theta];
fi;
od;
od;
if numelems(good_d_found) = 0 then return FAIL; fi;
#now convert the array to Matrix, and sort on d, so
#we start will the smallest degree d, which is the first column, as
#that will be most efficient.
#convert to set first, to remove any possible duplicat entries
#then convert to Matrix
convert(good_d_found,set);
B := convert(convert(%,list),Matrix);
B := B[sort(B[.., 1], 'output'= 'permutation')];
for N from 1 to LinearAlgebra:-RowDimension(B) do
r_solution := _self:-case_2_step_3(B[N,1], B[N,2]); #(d,theta)
if r_solution <> FAIL then
y_solution := _self:-build_y_solution_from_r_solution(r_solution);
return y_solution;
fi;
od;
return FAIL;
end proc;
#--------------------------------------------------------------------
# called from _self:-case_2_step_2()
#--------------------------------------------------------------------
local case_2_step_3::static:=proc(_self,d::nonnegative,theta,$)
local p, i;
local a::nothing;
local tmp;
local coeff_sol := [];
local eq;
local phi;
local sol_w;
local current_w;
local r:=_self:-r;
local w::nothing;
local x := _self:-x;
local sol;
p := x^d;
for i from d-1 by -1 to 0 do
p := p + a[i] * x^i;
od;
eq:= simplify(diff(p,x$3) + 3*theta*diff(p,x$2) +
(3*theta^2+3*diff(theta,x) -4*r) * diff(p,x)
+ ( diff(theta,x$2)+3*theta*diff(theta,x)+theta^3 -
4*r*theta - 2*diff(r,x))*p) = 0;
if d=0 then
if not evalb(eq) then return FAIL; fi;
else
#solve for polynomial coefficients
try
coeff_sol:= timelimit(30,solve(
identity(eq,x), [seq(a[i],i=0..d-1)]));
catch:
return FAIL;
end try;
if nops(coeff_sol) = 0 then return FAIL; fi;
tmp := map(evalb,coeff_sol[1]);
if has(tmp,true) then return FAIL; fi;
p := eval(p,coeff_sol[1]); #to force a[i] solutions to update
fi;
phi := theta + diff(p,x)/p;
eq := w^2 - phi*w + simplify(1/2*diff(phi,x)+1/2*phi^2-r) = 0;
try
sol_w := timelimit(30,solve(eq,[w])); #changed to []
sol_w := ListTools:-Flatten(sol_w);
catch:
return FAIL;
end try;
if nops(sol_w) = 0 then return FAIL; fi;
for current_w in sol_w do
current_w := radsimp(rhs(current_w));
#verify w before using it. Added 1/12/2022 4 PM
tmp := diff(current_w,x)+current_w^2;
if evalb( tmp= r) or is(tmp = r) or simplify(tmp-r)=0 then
try
sol := timelimit(30,int(current_w, x));
catch:
return FAIL;
end try;
if has(sol,int) then return FAIL; fi;
return simplify(exp(sol));
fi;
od;
end proc;
#--------------------------------------------------------------------
#
# C A S E T H R E E I M P L E M E N T A T I O N
#
# returns ode solution using case 3, or FAIL is no solution exist
#--------------------------------------------------------------------
local solve_case_3::static:=proc(_self,$)
local E_inf::set;
local gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry):={};
local sol;
#these are possible degress of w to try until one works or none works
#local w_degree::list := [4,6,12];
local w_degree::list := [4,6,12];
local n::posint;
for n in w_degree do
gamma_set, E_inf := _self:-case_3_step_1(n);
sol := _self:-case_3_step_2(gamma_set, E_inf, n );
if sol<>FAIL then return sol; fi;
od;
return FAIL;
end proc;
#--------------------------------------------------------------------
# First step in case 3
#--------------------------------------------------------------------
local case_3_step_1::static:=proc(_self,
n::posint,
$)::set(kovacic_class:-case_2_and_3_gamma_entry),set;
local current_pole;
local e::kovacic_class:-case_2_and_3_gamma_entry;
local E_inf::set;
local b,k;
local x := _self:-x;
local r := _self:-r;
#this contains all information generated for each pole of r
#this is what is called the set GAMMA in the paper and in the diagram of
#algorithm above.
local gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry):={};
for current_pole in _self:-poles_list do
e := Object(kovacic_class:-case_2_and_3_gamma_entry);
e:-pole_location := current_pole[1];
e:-pole_order := current_pole[2];
if e:-pole_order = 1 then
e:-Ec := {12};
gamma_set := { op(gamma_set), e };
elif e:-pole_order = 2 then
e:-b := _self:-b_partial_fraction(r,x,e:-pole_location,2);
e:-Ec := {seq( 6+12*k/n*sqrt(1+4*(e:-b)),k=-n/2..n/2,1)};
e:-Ec := select(z->type(z,integer),e:-Ec);
gamma_set := { op(gamma_set), e };
else
ERROR("Internal error. Case 3 can only have poles of order 1 or 2");
fi;
od;
#same formula, but different b
b := lcoeff(_self:-s) / lcoeff(_self:-t);
E_inf := {seq( 6+12*k/n*sqrt(1+4*b),k=-n/2..n/2,1)};
E_inf := select(z->type(z,integer),E_inf);
return gamma_set,E_inf;
end proc;
#--------------------------------------------------------------------
# Second step in case 3
#
# called from _self:-solve_case_3()
# This determines the set of d non-negative integers and
# corresponding w for each d.
#--------------------------------------------------------------------
local case_3_step_2::static:=proc(
_self,
gamma_set::set(kovacic_class:-case_2_and_3_gamma_entry),
E_inf::set,
n::posint,$)
local L::list := [];
local item,current_E_inf;
local d,ee,theta,N;
local x := _self:-x;
local r_solution,y_solution;
local S;
local current_iteration::integer;
local tmp;
#this will contains all data found for any nonnegative d. Each
#entry will be a list of this form
# [d, theta, S]
#so if we obtain say 3 values of d that are nonnegative
local good_d_found := Array(1..0);
local B::Matrix; #this will contain good_d_found as matrix
#DEBUG();
for item in gamma_set do
L := [ op(L), convert(item:-Ec,list) ];
od;
#now find all possible tuples
if nops(L)>1 then
L := kovacic_class:-cartProdSeq(op(L));
else
L := L[1];
fi;
current_iteration:=0;
for current_E_inf in E_inf do
for item in L do
current_iteration := current_iteration+1;
d := n/12*( current_E_inf - add(item));
if type(d,'integer') and d >= 0 then
theta :=0;
for N,ee in item do
theta := theta + ee/(x-gamma_set[N]:-pole_location);
od;
theta := n/12*theta;
theta := simplify(theta);
S := mul( (x-gamma_set[N]:-pole_location), N=1..nops(item));
tmp := simplify(S) assuming real;
if not has(tmp,csgn) and not has(tmp,signum) then
S:= tmp;
fi;
good_d_found ,= [ d, theta, S];
fi;
od;
od;
#now convert the array to Matrix, and sort on d, so we start
# with the smallest
#degree d, which is the first column, as that will be most efficient.
if numelems(good_d_found) = 0 then return FAIL; fi;
#convert to set first, to remove any possible duplicat entries
#then convert to Matrix
convert(good_d_found,set);
B := convert(convert(%,list),Matrix);
B := B[sort(B[.., 1], 'output'= 'permutation')];
for N from 1 to LinearAlgebra:-RowDimension(B) do
#(d,theta,S,n)
r_solution := _self:-case_3_step_3(B[N,1], B[N,2], B[N,3],n);
if r_solution <> FAIL then
_self:-n_used_for_case_3 := n;
y_solution := _self:-build_y_solution_from_r_solution(r_solution);
return y_solution;
fi;
od;
return FAIL;
end proc;
#--------------------------------------------------------------------
# Third and final step in case 3.
# called from _self:-case_3_step_2(). Returns solution or FAIL is no
# solution found.
#--------------------------------------------------------------------
local case_3_step_3::static:=proc(
_self,
d::nonnegative,
theta,
S,
n::posint,
$)
local p, i,P_minus_1;
local a::nothing;
local final_result,tmp;
local coeff_sol := [];
local sol_w;
local current_w;
local r := _self:-r;
local x := _self:-x;
local omega::symbol;
local P := Array(-1..n);
local result_of_simplify;
local result_of_is_check;
local omega_equation;
#this makes p(x). For example if d=3, then the result will be
# p(x) = x^3 + a(2) x^2 + a(1) x + a(0)
#where the number of unknowns to determine is always the same
#as the degree, in this case a(0),a(1),a(2)
p := x^d; #this will be 1 if degree is zero
for i from d-1 by -1 to 0 do
p := p + a[i] * x^i;
od;
#build the P_n polynomials based on p(x) above
P[n] := -p;
for i from n by -1 to 0 do
if i=n then
P[i-1] := -S * diff(P[i], x) - S*theta*P[i];
P[i-1] := simplify( P[i-1]);
else
P[i-1] := -S * diff(P[i], x) +
( (n-i)*diff(S,x) - S*theta)*P[i] - (n-i)*(i+1)*S^2*r*P[i+1];
P[i-1] := simplify( P[i-1]);
fi;
od;
#solve P[-1] for p(x) if needed.
P_minus_1 := expand(numer(radsimp(P[-1])));
if P_minus_1 <>0 and evalb(p<>1) then
if not hastype(P_minus_1,'indexed') then #it must be indexed now
ERROR("Internal error. Please report. This should not happen");
fi;
try
coeff_sol := timelimit(30,solve(identity(P_minus_1,x), [seq(a[i],i=0..d-1)]));
catch:
return FAIL;
end try;
if nops(coeff_sol)=0 then #unable to solve
return FAIL;
fi;
map(evalb,coeff_sol[1]); #check all solved for
if has(%,true) then return FAIL; fi;
fi;
#build the equation for omega
omega_equation := 0;
for i from 0 to n do
omega_equation := omega_equation + S^i*P[i]/(n-i)! * omega^i ;
od;
omega_equation := simplify(omega_equation);
if nops(coeff_sol)<>0 then
#to force a[i] solutions to update to solved coefficients
omega_equation := eval(omega_equation,coeff_sol[1]);
fi;
try
sol_w := timelimit(30, solve(omega_equation=0, [omega]));
sol_w := ListTools:-Flatten(sol_w);
catch:
return FAIL;
end try;
if nops(sol_w) = 0 then return FAIL; fi;
#go over each w solution and use one that works. verify before using
for current_w in sol_w do
current_w := rhs(current_w);
if not has(current_w, RootOf) then
try
current_w := timelimit(30,simplify(current_w));
catch:
NULL;
end try;
if is_w_verified(current_w,x,r) then
final_result := _self:-simplify_final_result(current_w,x);
if final_result<>FAIL then
return final_result;
fi;
fi;
else #no Rootof, try to resolve
try
current_w := timelimit(30,[allvalues(current_w)]);
catch:
return FAIL;
end try;
current_w := current_w[1]; #just use any root. Pick first
if not has(current_w, RootOf) then
try
current_w := timelimit(30,simplify(current_w));
catch:
NULL;
end try;
if is_w_verified(current_w,x,r) then
final_result := _self:-simplify_final_result(current_w,x);
if final_result<>FAIL then
return final_result;
fi;
fi;
fi;
fi;
od;
return FAIL;
end proc;
#--------------------------------------------------------------------
# Called from case 3, step 3 to verify w
#--------------------------------------------------------------------
local is_w_verified:=proc(current_w,x,r)::truefalse;
local tmp;
local result_of_is_check::truefalse;
local result_of_simplify;
tmp := diff(current_w,x)+current_w^2;
if evalb( tmp=r) then return true; fi;
try
result_of_simplify := timelimit(30,simplify(tmp-r));
if evalb(result_of_simplify=0) then return true; fi;
catch:
NULL;
end try;
try
result_of_is_check := timelimit(30,is(tmp = r));
if result_of_is_check then return true; fi;
catch:
NULL;
end try;
return false;
end proc;
#--------------------------------------------------------------------
# Called from case 3, step 3 to simplify final result
#--------------------------------------------------------------------
local simplify_final_result::static:=proc(_self,omega,x,$)
local final_result;
local integral_result;
try
integral_result := timelimit(30,int(omega, x));
catch:
return FAIL;
end try;
if has(integral_result,int) or has(integral_result,Int) then
return exp(Int(omega, x));
fi;
try
final_result := timelimit(30,simplify(exp(integral_result))) assuming real;
if has(final_result,signum) or has(final_result,csgn)
or has(final_result,abs) then
final_result := exp(integral_result);
fi;
catch:
final_result:= exp(integral_result);
end try;
return final_result;
end proc;
#--------------------------------------------------------------------
# External proc helper
# provided thanks to Joseph Riel
#--------------------------------------------------------------------
local cartProdSeq:= proc(L::seq(list))
local Seq::nothing,i::nothing,j;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
eval([subs(Seq= seq, foldl(Seq, [cat(i, 1..nargs)],
seq(cat(i,j)= L[j], j= nargs..1, -1)))])
end proc:
#--------------------------------------------------------------------
# Checks for special math function.
# This function was provided thanks to Carl Love.
# modified to allow some special functions.
#--------------------------------------------------------------------
local has_special_math_function:= subs(
_F= {(op@FunctionAdvisor)~(FunctionAdvisor("class_members", "quiet"), "quiet")[]}
minus ( {FunctionAdvisor("elementary", "quiet")[]}
union {erf,erfc, erfi, Im,Re,signum,max,argument} ),
proc(e::algebraic, x::{name, set(name)}:= {})
hastype(e, And(specfunc(_F), dependent(x)))
end proc
):
#--------------------------------------------------------------------
# called if step 3 is successful. Build y solution from r solution
#--------------------------------------------------------------------
local build_y_solution_from_r_solution::static:=proc(_self,r_solution,$)
local y1,y2;
local int_B_over_A;
local A :=_self:-A, B:=_self:-B;
local x :=_self:-x, y:= _self:-y;
local tmp;
if B = 0 then
y1 := r_solution;
int_B_over_A := 0;
else
try
int_B_over_A := timelimit(60,int(-B/A,x));
if has(int_B_over_A,int) or has(int_B_over_A,RootOf) then
int_B_over_A := Int(-B/A, x);
y1 :=r_solution*exp(1/2*int_B_over_A);
else
y1 := simplify(r_solution*exp(1/2*int_B_over_A));
if has(y1,signum) or has(y1,csgn) then
y1 := r_solution*exp(1/2*int_B_over_A);
fi;
fi;
catch:
int_B_over_A := Int(-B/A,x);
y1 := r_solution*exp(1/2*int_B_over_A);
end try;
fi;
try
y2 := timelimit(30, int( simplify(exp(int_B_over_A))/y1^2,x));
if kovacic_class:-has_special_math_function(y2,x) then
y2 := y1 * Int( exp(int_B_over_A)/y1^2,x)
else
y2 := y1 * y2;
fi;
tmp := simplify(y2);
if has(tmp,signum) or has(tmp,csgn) then
y2 := y1 * y2;
else
y2 := tmp;
fi;
catch:
y2 := y1 * Int( exp(int_B_over_A)/y1^2,x)
end try;
return y(x) = _C1*y1 + _C2*y2;
end proc;
#--------------------------------------------------------------------
# Helper private function called to obtain b from the partial fraction
# decomposition of r needed by the differenent case implementations
#--------------------------------------------------------------------
local b_partial_fraction::static:=proc(_self,
r,
x::symbol,
pole_location,
the_power::posint,$)
local r_partial_fraction,T::nothing;
r_partial_fraction := allvalues(convert(r,'fullparfrac',x));
#adding dummy T so that select below always works
r_partial_fraction := T+r_partial_fraction;
select(z->hastype(z,
anything/(anything*identical(x - pole_location)^the_power)),
r_partial_fraction);
return coeff(%,1/(x-pole_location)^the_power);
end proc;
#--------------------------------------------------------------------
# Helper private function called to return Laurent series coefficient
# of the 1/(x-c)^n term. where n=1,2,3,.... for the function f(x)
# expandid about pole of order m
# if c=infinity, then x is replaced by 1/y and expansion is around 0
# for infinity, n=0,1,2,...
#--------------------------------------------------------------------
local laurent_coeff::static:=proc(_self,
f,
x::symbol,
c,
m::integer,
n::integer,$)
local the_coeff,y::nothing,fy;
if c = infinity then
fy := eval(f,x=1/y);
if m-n>0 then
the_coeff := limit( diff( y^m * fy,y$(m-n)),y=0,right)/(m-n)!;
elif m-n=0 then
the_coeff := limit( y^m * fy,y=0,right)/(m-n)!;
else
the_coeff := 0;
fi;
else
if m-n>0 then
the_coeff := limit( diff( (x-c)^m * f,x$(m-n)),x=c,right)/(m-n)! ;
elif m-n=0 then
the_coeff := limit( (x-c)^m * f,x=c,right)/(m-n)!;
else
the_coeff := 0;
fi;
fi;
the_coeff := eval(the_coeff,[csgn=1,signum=1]);
return the_coeff;
end proc;
end module;