1 /* -*- mode:C++ ; compile-command: "g++-3.4 -I.. -g -c desolve.cc -DHAVE_CONFIG_H -DIN_GIAC" -*- */ 2 #include "giacPCH.h" 3 /* 4 * Copyright (C) 2000, 2014 B. Parisse, Institut Fourier, 38402 St Martin d'Heres 5 * 6 * This program is free software; you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation; either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. 18 */ 19 using namespace std; 20 #include <stdexcept> 21 #include <cmath> 22 #include "desolve.h" 23 #include "derive.h" 24 #include "intg.h" 25 #include "subst.h" 26 #include "usual.h" 27 #include "symbolic.h" 28 #include "unary.h" 29 #include "poly.h" 30 #include "sym2poly.h" // for equalposcomp 31 #include "tex.h" 32 #include "modpoly.h" 33 #include "series.h" 34 #include "solve.h" 35 #include "ifactor.h" 36 #include "prog.h" 37 #include "rpn.h" 38 #include "lin.h" 39 #include "intgab.h" 40 #include "giacintl.h" 41 42 #ifndef NO_NAMESPACE_GIAC 43 namespace giac { 44 #endif // ndef NO_NAMESPACE_GIAC 45 integrate_without_lnabs(const gen & e,const gen & x,GIAC_CONTEXT)46 gen integrate_without_lnabs(const gen & e,const gen & x,GIAC_CONTEXT){ 47 // workaround for desolve(diff(y)*sin(x)=y*ln(y),x,y); 48 // otherwise it returns ln(-1-cos(x)) 49 bool save_cv=complex_variables(contextptr); 50 complex_variables(false,contextptr); 51 gen res=integrate_gen(e,x,contextptr); 52 if (lop(res,at_abs).empty() && lop(res,at_floor).empty()){ 53 complex_variables(save_cv,contextptr); 54 return res; 55 } 56 bool save_do_lnabs=do_lnabs(contextptr); 57 do_lnabs(false,contextptr); 58 res=integrate_gen(e,x,contextptr); 59 do_lnabs(save_do_lnabs,contextptr); 60 complex_variables(save_cv,contextptr); 61 return res; 62 } 63 gen_t(const vecteur & v,GIAC_CONTEXT)64 gen gen_t(const vecteur & v,GIAC_CONTEXT){ 65 #ifdef GIAC_HAS_STO_38 66 identificateur id_t("t38_"); 67 #else 68 identificateur id_t(" t"); 69 #endif 70 gen tmp_t,t=t__IDNT; 71 t=t._IDNTptr->eval(1,tmp_t,contextptr); 72 if (t!=t__IDNT || equalposcomp(lidnt(v),t__IDNT)) 73 t=id_t; 74 return t; 75 } 76 laplace(const gen & f0,const gen & x,const gen & s,GIAC_CONTEXT)77 gen laplace(const gen & f0,const gen & x,const gen & s,GIAC_CONTEXT){ 78 if (x.type!=_IDNT) 79 return gensizeerr(contextptr); 80 if (f0.type==_VECT){ 81 vecteur v=*f0._VECTptr; 82 for (int i=0;i<int(v.size());++i){ 83 v[i]=laplace(v[i],x,s,contextptr); 84 } 85 return gen(v,f0.subtype); 86 } 87 gen t(s); 88 if (s==x){ 89 #ifdef GIAC_HAS_STO_38 90 t=identificateur("s38_"); 91 #else 92 t=identificateur(" t"); 93 #endif 94 } 95 // check for negative powers of x in f 96 gen f(f0); 97 vecteur v(1,x); 98 lvar(f,v); 99 fraction ff=sym2r(f,v,contextptr); 100 gen ffden=ff.den; 101 int n=0; 102 if (ffden.type==_POLY){ 103 polynome & ffdenp = *ffden._POLYptr; 104 if (!ffdenp.coord.empty() && (n=ffdenp.coord.back().index.front()) ){ 105 // multiply by (-1)^n*x^n, do laplace, then integrate n times answer 106 index_t idxt(v.size()); 107 idxt.front()=-n; 108 ff=fraction(ff.num,ffden._POLYptr->shift(idxt)); 109 f=r2sym(ff,v,contextptr); 110 if (n%2) 111 f=-f; 112 } 113 } 114 if (!assume_t_in_ab(t,plus_inf,plus_inf,true,true,contextptr)) 115 return gensizeerr(contextptr); 116 gen res=_integrate(makesequence(f*exp(-t*x,contextptr),x),contextptr); 117 if (lop(res,at_integrate).empty()) 118 res=-_limit(makesequence(res,x,0,1),contextptr); 119 else 120 res=undef; 121 if (is_undef(res)) 122 res=_integrate(makesequence(f*exp(-t*x,contextptr),x,0,plus_inf),contextptr); 123 for (int i=1;i<=n;++i){ 124 if (is_undef(res)) 125 return res; 126 res = _integrate(gen(makevecteur(res,t,0,t),_SEQ__VECT),contextptr); 127 res += _integrate(gen(makevecteur(f/pow(-x,i),x,0,plus_inf),_SEQ__VECT),contextptr); 128 } 129 purgenoassume(t,contextptr); 130 if (s==x) 131 res=subst(res,t,x,false,contextptr); 132 return ratnormal(res,contextptr); 133 /* 134 gen remains,res=integrate(f*exp(-t*x,contextptr),*x._IDNTptr,remains,contextptr); 135 res=subst(-res,x,zero,false,contextptr); 136 if (s==x) 137 res=subst(res,t,x,false,contextptr); 138 if (!is_zero(remains)) 139 res = res +symbolic(at_integrate,gen(makevecteur(remains,x,0,plus_inf),_SEQ__VECT)); 140 return res; 141 */ 142 } 143 _laplace_(const gen & args,GIAC_CONTEXT)144 static gen _laplace_(const gen & args,GIAC_CONTEXT){ 145 if (args.type!=_VECT) 146 return laplace(args,vx_var,vx_var,contextptr); 147 vecteur & v=*args._VECTptr; 148 int s=int(v.size()); 149 if (s==2) 150 return laplace( v[0],v[1],v[1],contextptr); 151 if (s!=3) 152 return gensizeerr(contextptr); 153 return laplace( v[0],v[1],v[2],contextptr); 154 } 155 // "unary" version _laplace(const gen & args,GIAC_CONTEXT)156 gen _laplace(const gen & args,GIAC_CONTEXT){ 157 if ( args.type==_STRNG && args.subtype==-1) return args; 158 bool b=approx_mode(contextptr); 159 approx_mode(false,contextptr); 160 #ifndef NSPIRE 161 my_ostream * ptr=logptr(contextptr); 162 logptr(0,contextptr); 163 gen res=_laplace_(args,contextptr); 164 logptr(ptr,contextptr); 165 #else 166 gen res=_laplace_(exact(args,contextptr),contextptr); 167 #endif 168 approx_mode(b,contextptr); 169 if (b || has_num_coeff(args)) 170 res=simplifier(evalf(res,1,contextptr),contextptr); 171 return res; 172 } 173 static const char _laplace_s []="laplace"; 174 static define_unary_function_eval (__laplace,&_laplace,_laplace_s); 175 define_unary_function_ptr5( at_laplace ,alias_at_laplace,&__laplace,0,true); 176 cstcoeff(const polynome & p)177 polynome cstcoeff(const polynome & p){ 178 vector< monomial<gen> >::const_iterator it=p.coord.begin(),itend=p.coord.end(); 179 for (;it!=itend;++it){ 180 if (it->index.front()==0) 181 break; 182 } 183 return polynome(p.dim,vector< monomial<gen> >(it,itend)); 184 } 185 186 // reduction of a fraction with multiple poles to single poles by integration 187 // by part, use the relation 188 // ilaplace(P'/P^(k+1))=laplacevar/k*ilaplace(1/P^k) laplace_reduce_pf(const pf<gen> & p_cst,tensor<gen> & laplacevar)189 pf<gen> laplace_reduce_pf(const pf<gen> & p_cst, tensor<gen> & laplacevar ){ 190 pf<gen> p(p_cst); 191 assert(p.mult>0); 192 if (p.mult==1) 193 return p_cst; 194 tensor<gen> fprime=p.fact.derivative(); 195 tensor<gen> d(fprime.dim),C(fprime.dim),u(fprime.dim),v(fprime.dim); 196 egcdpsr(p.fact,fprime,u,v,d); // f*u+f'*v=d 197 tensor<gen> usave(u),vsave(v); 198 // int initial_mult=p.mult-1; 199 while (p.mult>1){ 200 egcdtoabcuv(p.fact,fprime,p.num,u,v,d,C); 201 p.mult--; 202 p.den=(p.den/p.fact)*C*gen(p.mult); 203 p.num=u*gen(p.mult)+v.derivative()+v*laplacevar; 204 if ( (p.mult % 5)==1) // simplify from time to time 205 TsimplifybyTlgcd(p.num,p.den); 206 if (p.mult==1) 207 break; 208 u=usave; 209 v=vsave; 210 } 211 return pf<gen>(p); 212 } 213 pf_ilaplace(const gen & e0,const gen & x,gen & remains,GIAC_CONTEXT)214 static gen pf_ilaplace(const gen & e0,const gen & x, gen & remains,GIAC_CONTEXT){ 215 vecteur vexp; 216 gen res; 217 lin(e0,vexp,contextptr); // vexp = coeff, arg of exponential 218 const_iterateur it=vexp.begin(),itend=vexp.end(); 219 remains=0; 220 for (;it!=itend;){ 221 gen coeff=*it; 222 ++it; 223 gen axb=*it,expa,expb; 224 ++it; 225 gen e=coeff*exp(axb,contextptr); 226 if (!is_linear_wrt(axb,x,expa,expb,contextptr)){ 227 remains += e; 228 continue; 229 } 230 if (is_strictly_positive(expa,contextptr)) 231 *logptr(contextptr) << gettext("Warning, exponential x coeff is positive ") << expa << '\n'; 232 vecteur varx(lvarx(coeff,x)); 233 int varxs=int(varx.size()); 234 if (!varxs){ // Dirac function 235 res += coeff*exp(expb,contextptr)*symbolic(at_Dirac,laplace_var+expa); 236 continue; 237 } 238 if ( (varxs>1) || (varx.front()!=x) ) { 239 remains += e; 240 continue; 241 } 242 vecteur l; 243 l.push_back(x); // insure x is the main var 244 l.push_back(laplace_var); // s var as second var 245 l=vecteur(1,l); 246 alg_lvar(makevecteur(coeff,axb),l); 247 gen glap=e2r(laplace_var,l,contextptr); 248 if (glap.type!=_POLY) return gensizeerr(gettext("desolve.cc/pf_ilaplace")); 249 int s=int(l.front()._VECTptr->size()); 250 if (!s){ 251 l.erase(l.begin()); 252 s=int(l.front()._VECTptr->size()); 253 } 254 gen r=e2r(coeff,l,contextptr); 255 gen r_num,r_den; 256 fxnd(r,r_num,r_den); 257 if (r_num.type==_EXT){ 258 remains += e; 259 continue; 260 } 261 if (r_den.type!=_POLY){ 262 remains += e; 263 continue; 264 } 265 polynome den(*r_den._POLYptr),num(s); 266 if (r_num.type==_POLY) 267 num=*r_num._POLYptr; 268 else 269 num=polynome(r_num,s); 270 polynome p_content(lgcd(den)); 271 factorization vden(sqff(den/p_content)); // first square-free factorization 272 vector< pf<gen> > pfde_VECT; 273 polynome ipnum(s),ipden(s),temp(s),tmp(s); 274 partfrac(num,den,vden,pfde_VECT,ipnum,ipden); 275 vector< pf<gen> >::iterator it=pfde_VECT.begin(); 276 vector< pf<gen> >::const_iterator itend=pfde_VECT.end(); 277 vector< pf<gen> > rest,finalde_VECT; 278 for (;it!=itend;++it){ 279 pf<gen> single(laplace_reduce_pf(*it,*glap._POLYptr)); 280 gen extra_div=1; 281 factor(single.den,p_content,vden,false,withsqrt(contextptr),complex_mode(contextptr),1,extra_div); 282 partfrac(single.num,single.den,vden,finalde_VECT,temp,tmp); 283 } 284 it=finalde_VECT.begin(); 285 itend=finalde_VECT.end(); 286 gen lnpart(0),deuxaxplusb,sqrtdelta,exppart; 287 polynome a(s),b(s),c(s); 288 polynome d(s),E(s),lnpartden(s); 289 polynome delta(s),atannum(s),alpha(s); 290 vecteur lprime(l); 291 if (lprime.front().type!=_VECT) return gensizeerr(gettext("desolve.cc/pf_ilaplace")); 292 lprime.front()=cdr_VECT(*(lprime.front()._VECTptr)); 293 bool uselog; 294 for (;it!=itend;++it){ 295 int deg=it->fact.lexsorted_degree(); 296 switch (deg) { 297 case 1: // 1st order 298 findde(it->den,a,b); 299 lnpart=lnpart+rdiv(r2e(it->num,l,contextptr),r2e(firstcoeff(a),lprime,contextptr),contextptr)*exp(r2e(rdiv(-b,a,contextptr),lprime,contextptr)*laplace_var,contextptr); 300 break; 301 case 2: // 2nd order 302 findabcdelta(it->fact,a,b,c,delta); 303 exppart=exp(r2e(rdiv(-b,gen(2)*a,contextptr),lprime,contextptr)*laplace_var,contextptr); 304 uselog=is_positive(delta); 305 alpha=(it->den/it->fact).trunc1()*a; 306 findde(it->num,d,E); 307 atannum=a*E*gen(2)-b*d; 308 // cos part d/alpha*ln(fact) 309 lnpartden=alpha; 310 simplify(d,lnpartden); 311 if (uselog){ 312 sqrtdelta=normal(sqrt(r2e(delta,lprime,contextptr),contextptr),contextptr); 313 gen racine=ratnormal(sqrtdelta/gen(2)/r2e(a,lprime,contextptr),contextptr); 314 lnpart=lnpart+rdiv(r2e(d,lprime,contextptr),r2e(lnpartden,lprime,contextptr),contextptr)*cosh(racine*laplace_var,contextptr)*exppart; 315 gen aa=ratnormal(r2e(atannum,lprime,contextptr)/r2e(alpha,lprime,contextptr)/sqrtdelta,contextptr); 316 lnpart=lnpart+aa*sinh(racine*laplace_var,contextptr)*exppart; 317 } 318 else { 319 sqrtdelta=normal(sqrt(r2e(-delta,lprime,contextptr),contextptr),contextptr); 320 gen racine=ratnormal(sqrtdelta/gen(2)/r2e(a,lprime,contextptr),contextptr); 321 lnpart=lnpart+rdiv(r2e(d,lprime,contextptr),r2e(lnpartden,lprime,contextptr),contextptr)*cos(racine*laplace_var,contextptr)*exppart; 322 gen aa=ratnormal(r2e(atannum,lprime,contextptr)/r2e(alpha,lprime,contextptr)/sqrtdelta,contextptr); 323 lnpart=lnpart+aa*sin(racine*laplace_var,contextptr)*exppart; 324 } 325 break; 326 default: 327 rest.push_back(pf<gen>(it->num,it->den,it->fact,1)); 328 break ; 329 } 330 } 331 vecteur ipnumv=polynome2poly1(ipnum,1); 332 gen deno=r2e(ipden,l,contextptr); 333 int nums=int(ipnumv.size()); 334 for (int i=0;i<nums;++i){ 335 gen tmp = rdiv(r2e(ipnumv[i],lprime,contextptr),deno,contextptr); 336 tmp = tmp*symbolic(at_Dirac,(i==nums-1)?laplace_var:gen(makevecteur(laplace_var,nums-1-i),_SEQ__VECT)); 337 res += tmp; 338 } 339 remains += r2sym(rest,l,contextptr)*exp(axb,contextptr); 340 if (is_zero(expa)) 341 res += lnpart*exp(expb,contextptr); 342 else 343 res += quotesubst(lnpart,laplace_var,laplace_var+expa,contextptr)*exp(expb,contextptr)*_Heaviside(laplace_var+expa,contextptr); 344 } 345 return res; 346 } 347 ilaplace(const gen & f,const gen & x,const gen & s,GIAC_CONTEXT)348 gen ilaplace(const gen & f,const gen & x,const gen & s,GIAC_CONTEXT){ 349 if (x.type!=_IDNT) 350 return gensizeerr(contextptr); 351 if (has_num_coeff(f)) 352 return ilaplace(exact(f,contextptr),x,s,contextptr); 353 gen remains,res=linear_apply(f,x,remains,contextptr,pf_ilaplace); 354 res=subst(res,laplace_var,s,false,contextptr); 355 if (!is_zero(remains)) 356 res=res+symbolic(at_ilaplace,makevecteur(remains,x,s)); 357 return res; 358 } 359 // "unary" version _ilaplace(const gen & args,GIAC_CONTEXT)360 gen _ilaplace(const gen & args,GIAC_CONTEXT){ 361 if ( args.type==_STRNG && args.subtype==-1) return args; 362 if (args.type!=_VECT) 363 return ilaplace(args,vx_var,vx_var,contextptr); 364 vecteur & v=*args._VECTptr; 365 int s=int(v.size()); 366 if (s==2) 367 return ilaplace( v[0],v[1],v[1],contextptr); 368 if (s!=3) 369 return gensizeerr(contextptr); 370 return ilaplace( v[0],v[1],v[2],contextptr); 371 } 372 static const char _ilaplace_s []="ilaplace"; 373 static define_unary_function_eval (__ilaplace,&_ilaplace,_ilaplace_s); 374 define_unary_function_ptr5( at_ilaplace ,alias_at_ilaplace,&__ilaplace,0,true); 375 376 static const char _invlaplace_s []="invlaplace"; 377 static define_unary_function_eval (__invlaplace,&_ilaplace,_invlaplace_s); 378 define_unary_function_ptr5( at_invlaplace ,alias_at_invlaplace,&__invlaplace,0,true); 379 unable_to_solve_diffeq()380 static gen unable_to_solve_diffeq(){ 381 return gensizeerr(gettext("Unable to solve differential equation")); 382 } 383 diffeq_constante(int i,GIAC_CONTEXT)384 gen diffeq_constante(int i,GIAC_CONTEXT){ 385 #if 0 // def NSPIRE 386 if (i<5){ 387 const char * tab[]={"o","p","q","r","s"}; 388 return gen(tab[i],contextptr); 389 } 390 #endif 391 #ifdef GIAC_HAS_STO_38 392 string s("G_"+print_INT_(i)); 393 #else 394 string s("c_"+print_INT_(i)); 395 #endif 396 return gen(s,contextptr); 397 } 398 399 // return -1 if f does not depend on y diffeq_order(const gen & f,const gen & y)400 static int diffeq_order(const gen & f,const gen & y){ 401 vecteur ydepend(rlvarx(f,y)); 402 const_iterateur it=ydepend.begin(),itend=ydepend.end(); 403 // since we did a recursive lvar we dismiss all variables except 404 // if they begin with derive 405 int n=-1; 406 for (;it!=itend;++it){ 407 if (*it==y) 408 n=giacmax(n,0); 409 if ( (it->type==_SYMB) && (it->_SYMBptr->sommet==at_derive) ){ 410 gen & g=it->_SYMBptr->feuille; 411 int m=-1,nder=1; 412 if ( (g.type==_VECT) && (!g._VECTptr->empty()) ){ 413 m=diffeq_order(g._VECTptr->front(),y); 414 if (g._VECTptr->size()==3){ 415 gen & gg=g._VECTptr->back(); 416 if (gg.type==_INT_) 417 nder=gg.val; 418 } 419 } 420 else 421 m=diffeq_order(g,y); 422 if (m>=0) 423 n=giacmax(n,m+nder); 424 } 425 } 426 return n; 427 } 428 429 // true if f is a linear differential equation 430 // & returns the coefficient in v in descending order 431 // v has size order+2 with last term=cst coeff of the diff equation is_linear_diffeq(const gen & f_orig,const gen & x,const gen & y,int order,vecteur & v,int step_info,GIAC_CONTEXT)432 static bool is_linear_diffeq(const gen & f_orig,const gen & x,const gen & y,int order,vecteur & v,int step_info,GIAC_CONTEXT){ 433 v.clear(); 434 gen f(f_orig),a,b,cur_y(y); 435 gen t=gen_t(makevecteur(x,y,f_orig),contextptr); 436 for (int i=0;i<=order;++i){ 437 gen ftmp(quotesubst(f,cur_y,t,contextptr)); 438 if (!is_linear_wrt(eval(ftmp,eval_level(contextptr),contextptr),t,a,b,contextptr)) 439 return false; 440 if (!rlvarx(a,y).empty()) 441 return false; 442 if (!i) 443 v.push_back(b); 444 v.push_back(a); 445 cur_y=symb_derive(y,x,i+1); 446 } 447 reverse(v.begin(),v.end()); 448 if (step_info && v.size()>3) 449 gprintf("Linear differential equation of coefficients %gen\nsecond member %gen",makevecteur(vecteur(v.begin(),v.end()-1),-v.back()),step_info,contextptr); 450 return true; 451 } 452 find_n_derivatives_function(const gen & f,const gen & x,int & nder,gen & fonction)453 static bool find_n_derivatives_function(const gen & f,const gen & x,int & nder,gen & fonction){ 454 if ( (f.type!=_SYMB) || (f._SYMBptr->sommet!=at_derive) ){ 455 nder=0; 456 fonction=f; 457 return true; 458 } 459 if (f._SYMBptr->feuille.type!=_VECT){ 460 if (!find_n_derivatives_function(f._SYMBptr->feuille,x,nder,fonction)) 461 return false; 462 ++nder; 463 return true; 464 } 465 vecteur & v=*f._SYMBptr->feuille._VECTptr; 466 if ( (v.size()>1) && (v[1]!=x) ) 467 return false; // setsizeerr(contextptr); 468 if (!find_n_derivatives_function(v[0],x,nder,fonction)) 469 return false; 470 if ( (v.size()==3) && (v[2].type==_INT_) ) 471 nder += v[2].val; 472 else 473 nder += 1; 474 return true; 475 } 476 function_of(const gen & y_orig,const gen & x_orig)477 static gen function_of(const gen & y_orig,const gen & x_orig){ 478 if ( (y_orig.type!=_SYMB) || (y_orig._SYMBptr->sommet!=at_of) ) 479 return gensizeerr(gettext("function_of")); 480 vecteur & v =*y_orig._SYMBptr->feuille._VECTptr; 481 if ( (v[1]!=x_orig) || (v[0].type!=_IDNT) ) 482 return gensizeerr(gettext("function_of")); 483 return v[0]; 484 } 485 in_desolve_with_conditions(const vecteur & v_,const gen & x,const gen & y,const gen & solution_generale,const vecteur & parameters,const gen & f,int step_info,GIAC_CONTEXT)486 static gen in_desolve_with_conditions(const vecteur & v_,const gen & x,const gen & y,const gen & solution_generale,const vecteur & parameters,const gen & f,int step_info,GIAC_CONTEXT){ 487 gen yy(y); 488 vecteur v(v_); 489 if (yy.type!=_IDNT) 490 yy=function_of(y,x); 491 if (is_undef(yy)) 492 return yy; 493 // special handling for systems 494 if (solution_generale.type==_VECT && v.size()==2){ 495 gen init=v[1],point=0; 496 if (init.is_symb_of_sommet(at_equal) && init._SYMBptr->feuille.type==_VECT&& init._SYMBptr->feuille._VECTptr->size()>=2){ 497 point=(*init._SYMBptr->feuille._VECTptr)[0]; 498 init=(*init._SYMBptr->feuille._VECTptr)[1]; 499 if (!point.is_symb_of_sommet(at_of) || point._SYMBptr->feuille.type!=_VECT || point._SYMBptr->feuille._VECTptr->size()<2 || point._SYMBptr->feuille._VECTptr->front()!=y) 500 return gensizeerr("Bad initial condition"); 501 point=(*point._SYMBptr->feuille._VECTptr)[1]; 502 } 503 gen systeme=subst(solution_generale,x,point,false,contextptr)-init; 504 gen s=_solve(makesequence(systeme,parameters),contextptr); 505 if (s.type!=_VECT) 506 return gensizeerr("Bad initial condition"); 507 vecteur res; 508 for (unsigned i=0;i<s._VECTptr->size();++i){ 509 gen tmp=subst(solution_generale,parameters,s[i],false,contextptr); 510 tmp=ratnormal(tmp,contextptr); 511 res.push_back(tmp); 512 } 513 return res; 514 } 515 if (solution_generale.type==_VECT) 516 *logptr(contextptr) << gettext("Boundary conditions for parametric curve not implemented") << '\n'; 517 // solve boundary conditions 518 iterateur jt=v.begin()+1,jtend=v.end(); 519 for (unsigned ndiff=0;jt!=jtend;++ndiff,++jt){ 520 if (jt->type==_VECT && jt->_VECTptr->size()==2){ 521 if (ndiff) 522 *jt=symbolic(at_of,makesequence(symbolic(at_derive,makesequence(y,x,int(ndiff))),jt->_VECTptr->front()))-jt->_VECTptr->back(); 523 else 524 *jt=symbolic(at_of,makesequence(y,jt->_VECTptr->front()))-jt->_VECTptr->back(); 525 } 526 } 527 const_iterateur it=v.begin()+1,itend=v.end(); 528 vecteur conditions(remove_equal(it,itend)); 529 if (conditions.empty()) 530 return solution_generale; 531 // conditions must be in terms of y(value) or derivatives 532 vecteur condvar(rlvarx(conditions,yy)); 533 vecteur yvar; // will contain triplet (var,n,x) n=nth derivative, x point 534 it=condvar.begin(),itend=condvar.end(); 535 int maxnder=0; 536 for (;it!=itend;++it){ 537 if ( (it->type!=_SYMB) || (it->_SYMBptr->sommet!=at_of) ) 538 continue; 539 vecteur & w=*it->_SYMBptr->feuille._VECTptr; 540 int nder; 541 gen fonction; 542 if (!find_n_derivatives_function(w[0],x,nder,fonction)) 543 return gensizeerr(contextptr); 544 if (fonction==y){ 545 if ( (w[1].type==_VECT) && (!w[1]._VECTptr->empty())) 546 yvar.push_back(makevecteur(*it,nder,w[1]._VECTptr->front())); 547 else 548 yvar.push_back(makevecteur(*it,nder,w[1])); 549 } 550 if (nder>maxnder) 551 maxnder=nder; 552 } 553 // compute all derivatives of the general solution 554 vecteur derivatives(1,solution_generale); 555 gen current=solution_generale; 556 for (int i=1;i<=maxnder;++i){ 557 current=derive(current,x,contextptr); 558 derivatives.push_back(current); 559 } 560 // evaluate at points of yvar making substition vectors 561 it=yvar.begin(),itend=yvar.end(); 562 vecteur substin,substout; 563 for (;it!=itend;++it){ 564 vecteur & w=*it->_VECTptr; 565 substin.push_back(w[0]); 566 substout.push_back(subst(derivatives[w[1].val],x,w[2],false,contextptr)); 567 } 568 // replace in conditions 569 conditions=*eval(subst(conditions,substin,substout,false,contextptr),eval_level(contextptr),contextptr)._VECTptr; 570 // solve system over _c0..._cn-1 571 int save_xcas_mode=xcas_mode(contextptr); 572 xcas_mode(contextptr)=0; 573 int save_calc_mode=calc_mode(contextptr); 574 calc_mode(contextptr)=0; 575 vecteur parameters_solutions=*_solve(gen(makevecteur(conditions,parameters),_SEQ__VECT),contextptr)._VECTptr; 576 if (step_info) 577 gprintf("General solution %gen\nSolving initial conditions\n%gen\nunknowns %gen\nSolutions %gen",makevecteur(solution_generale,conditions,parameters,parameters_solutions),step_info,contextptr); 578 xcas_mode(contextptr)=save_xcas_mode; 579 calc_mode(contextptr)=save_calc_mode; 580 // replace _c0..._cn-1 in solution_generale 581 it=parameters_solutions.begin(),itend=parameters_solutions.end(); 582 vecteur res; 583 for (;it!=itend;++it){ 584 gen solgen=eval(subst(solution_generale,parameters,*it,false,contextptr),eval_level(contextptr),contextptr); 585 // check if f is valid at points where conditions hold (3rd column of yvar) 586 gen solgenchk=eval(subst(f,y,solgen,false,contextptr),1,contextptr); 587 bool ok=true; 588 for (unsigned i=0;i<yvar.size();++i){ 589 gen tmp=subst(solgenchk,x,yvar[i][2],false,contextptr); 590 if (lidnt(tmp).empty() && !is_zero(simplify(tmp,contextptr))){ 591 ok=false; 592 break; 593 } 594 } 595 if (ok) 596 res.push_back(solgen); 597 } 598 if (res.size()==1) 599 return res.front(); 600 return res; 601 } 602 desolve_with_conditions(const vecteur & v,const gen & x,const gen & y,gen & f,int step_info,GIAC_CONTEXT)603 static gen desolve_with_conditions(const vecteur & v,const gen & x,const gen & y,gen & f,int step_info,GIAC_CONTEXT){ 604 if (v.empty()) 605 return gensizeerr(contextptr); 606 int ordre; bool num=false; 607 vecteur parameters; 608 gen solution_generale(desolve_f(v.front(),x,y,ordre,parameters,f,step_info,num,contextptr)); 609 if (solution_generale.type!=_VECT) { 610 gen res= in_desolve_with_conditions(v,x,y,solution_generale,parameters,f,step_info,contextptr); 611 return num?evalf(res,1,contextptr):res; 612 } 613 solution_generale.subtype=0; // otherwise desolve([y'=[[1,2],[2,1]]*y+[x,x+1],y(0)=[1,2]]) fails on the Prime (?) 614 if (parameters.empty()) 615 return num?evalf(solution_generale,1,contextptr):solution_generale; 616 iterateur it=solution_generale._VECTptr->begin(),itend=solution_generale._VECTptr->end(); 617 vecteur res; 618 res.reserve(itend-it); 619 for (;it!=itend;++it){ 620 if (it->type==_VECT) it->subtype=0; 621 gen tmp=in_desolve_with_conditions(v,x,y,*it,parameters,f,step_info,contextptr); 622 if (is_undef(tmp)) 623 return tmp; 624 if (tmp.type==_VECT) 625 res=mergevecteur(res,*tmp._VECTptr); 626 else 627 res.push_back(tmp); 628 } 629 return num?evalf(res,1,contextptr):res; 630 } 631 desolve_with_conditions(const vecteur & v,const gen & x,const gen & y,gen & f,GIAC_CONTEXT)632 static gen desolve_with_conditions(const vecteur & v,const gen & x,const gen & y,gen & f,GIAC_CONTEXT){ 633 int st=step_infolevel(contextptr); 634 step_infolevel(0,contextptr); 635 gen res=desolve_with_conditions(v,x,y,f,st,contextptr); 636 step_infolevel(st,contextptr); 637 return res; 638 } 639 640 // f must be a vector obtained using factors 641 // x, y are 2 idnt 642 // xfact and yfact should be initialized to 1 643 // return true if f=xfact*yfact where xfact depends on x and yfact on y only separate_variables(const gen & f,const gen & x,const gen & y,gen & xfact,gen & yfact,int step_info,GIAC_CONTEXT)644 bool separate_variables(const gen & f,const gen & x,const gen & y,gen & xfact,gen & yfact,int step_info,GIAC_CONTEXT){ 645 const_iterateur jt=f._VECTptr->begin(),jtend=f._VECTptr->end(); 646 for (;jt!=jtend;jt+=2){ 647 vecteur tmp(*_lname(*jt,contextptr)._VECTptr); 648 if (equalposcomp(tmp,y)){ 649 if (equalposcomp(tmp,x)) 650 return false; 651 yfact=yfact*pow(*jt,*(jt+1),contextptr); 652 } 653 else 654 xfact=xfact*pow(*jt,*(jt+1),contextptr); 655 } 656 if (step_info) 657 gprintf("Separable variables d%gen/%gen=%gen*d%gen",makevecteur(y,yfact,xfact,x),step_info,contextptr); 658 return true; 659 } 660 separate_variables(const gen & f,const gen & x,const gen & y,gen & xfact,gen & yfact,GIAC_CONTEXT)661 bool separate_variables(const gen & f,const gen & x,const gen & y,gen & xfact,gen & yfact,GIAC_CONTEXT){ 662 return separate_variables(f,x,y,xfact,yfact,step_infolevel(contextptr),contextptr); 663 } 664 ggb_varxy(const gen & f_orig,gen & vx,gen & vy,GIAC_CONTEXT)665 void ggb_varxy(const gen & f_orig,gen & vx,gen & vy,GIAC_CONTEXT){ 666 vecteur lv=lidnt(f_orig); 667 vx=vx_var; 668 vy=y__IDNT_e; 669 #if 0 670 if (calc_mode(contextptr)==1){ 671 vx=gen("ggbtmpvarx",contextptr); 672 vy=gen("ggbtmpvary",contextptr); 673 } 674 #endif 675 for (unsigned i=0;i<lv.size();++i){ 676 string s=lv[i].print(contextptr); 677 char c=s[s.size()-1]; 678 if (c=='x') 679 vx=lv[i]; 680 if (c=='y') 681 vy=lv[i]; 682 } 683 } 684 desolve_cleanup(const gen & i,const gen & x,GIAC_CONTEXT)685 static gen desolve_cleanup(const gen & i,const gen & x,GIAC_CONTEXT){ 686 if (i.is_symb_of_sommet(at_prod)){ 687 gen f=i._SYMBptr->feuille; 688 if (f.type==_VECT){ 689 vecteur w; 690 for (int j=0;j<f._VECTptr->size();++j){ 691 gen tmp=desolve_cleanup((*f._VECTptr)[j],x,contextptr); 692 if (!is_one(tmp)) 693 w.push_back(tmp); 694 } 695 return _prod(w,contextptr); 696 } 697 } 698 if (i.is_symb_of_sommet(at_abs) || i.is_symb_of_sommet(at_neg)) 699 return desolve_cleanup(i._SYMBptr->feuille,x,contextptr); 700 if (is_zero(derive(i,x,contextptr))) 701 return 1; 702 return i; 703 } 704 705 // solve linear diff eq of order 1 a*y'+b*y+c=0 desolve_lin1(const gen & a,const gen & b,const gen & c,const gen & x,vecteur & parameters,int step_info,GIAC_CONTEXT)706 static gen desolve_lin1(const gen &a,const gen &b,const gen & c,const gen & x,vecteur & parameters,int step_info,GIAC_CONTEXT){ 707 if (step_info) 708 gprintf("Linear differential equation of order 1 a*y'+b*y+c=0\na=%gen, b=%gen, c=%gen",makevecteur(a,b,c),step_info,contextptr); 709 if (a.type==_VECT){ 710 // y'+inv(a)*b(x)*y+inv(a)*c(x)=0 711 // take laplace transform 712 // p*Y-Y(0)+bsura*Y+csura=0 713 // (p+bsura)*Y=Y(0)-csura 714 int n=int(a._VECTptr->size()); 715 if (!ckmatrix(a) || !ckmatrix(b)) 716 return gensizeerr(contextptr); 717 gen inva=inv(a,contextptr); 718 gen bsura=inva*b,csura,cl; 719 if (!is_zero(derive(bsura,x,contextptr))) 720 return gensizeerr("Non constant linear differential system"); 721 if (c.type==_VECT){ 722 vecteur & cv=*c._VECTptr; 723 for (unsigned i=0;i<cv.size();++i){ 724 if (cv[i].type==_VECT && cv[i]._VECTptr->size()==1) 725 cv[i]=cv[i]._VECTptr->front(); 726 } 727 csura=inva*c; 728 cl=_laplace(makesequence(csura,x,x),contextptr); 729 } 730 else { 731 if (!is_zero(c)) 732 return gensizeerr("Invalid second member"); 733 cl=vecteur(n); 734 } 735 if (cl.type!=_VECT || int(cl._VECTptr->size())!=n) 736 return gensizeerr("Invalid second member"); 737 for (int i=0;i<n;++i){ 738 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 739 (*cl._VECTptr)[i] = parameters.back()- (*cl._VECTptr)[i]; 740 } 741 cl=inv(bsura+x,contextptr)*cl; 742 cl=ilaplace(cl,x,x,contextptr); 743 return vecteur(1,ratnormal(cl,contextptr)); 744 } 745 gen i=integrate_without_lnabs(inv(a,contextptr)*b,x,contextptr); 746 i=normal(lnexpand(i,contextptr),contextptr); 747 i=exp(i,contextptr); 748 if (step_info) 749 gprintf("Homogeneous solution C/%gen",makevecteur(i),step_info,contextptr); 750 i=expexpand(i,contextptr); 751 i=simplify(i,contextptr); 752 // cleanup general solution: remove cst factors and absolute values 753 i=desolve_cleanup(i,x,contextptr); 754 gen C=integrate_without_lnabs(ratnormal(rdiv(-c,a,contextptr)*i,contextptr),x,contextptr); 755 if (step_info && C!=0) 756 gprintf("Particuliar solution %gen",makevecteur(C),step_info,contextptr); 757 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 758 gen res=ratnormal(_lin((C+parameters.back())/i,contextptr),contextptr); 759 if (step_info) 760 gprintf("General solution %gen",makevecteur(res),step_info,contextptr); 761 return res; 762 } 763 desolve_linn(const gen & x,const gen & y,const gen & t,int n,vecteur & v,vecteur & parameters,gen & result,int step_info,GIAC_CONTEXT)764 bool desolve_linn(const gen & x,const gen & y,const gen & t,int n,vecteur & v,vecteur & parameters,gen & result,int step_info,GIAC_CONTEXT){ 765 // 1st order 766 if (n==1){ // a(x)*y'+b(x)*y+c(x)=0 767 // y'/y=-b/a -> y=C(x)exp(-int(b/a)) and a(x)*C'*exp()+c(x)=0 768 gen & a=v[0]; 769 gen & b=v[1]; 770 gen & c=v[2]; 771 if (ckmatrix(a)){ 772 if (c.type!=_VECT && is_zero(c)) 773 c=c*a; 774 c=_tran(c,contextptr)[int(a._VECTptr->size())-1]; 775 } 776 result=desolve_lin1(a,b,c,x,parameters,step_info,contextptr); 777 return true; 778 } 779 // cst coeff? 780 gen cst=v.back(); 781 v.pop_back(); 782 if (derive(v,x,contextptr)==vecteur(n+1,zero)){ 783 if (step_info) 784 gprintf("Linear differential equation with constant coefficients\nOrder %gen, coefficients %gen",makevecteur(n,v),step_info,contextptr); 785 // Yes! 786 // simpler general solution for small order generic lin diffeq with cst coeff/squarefree case 787 if (n<=3){ 788 vecteur rac=solve(horner(v,x,contextptr),x,1,contextptr); 789 comprim(rac); 790 if (n==2 && rac.size()==1){ 791 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 792 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 793 gen sol = exp(rac.front()*x,contextptr)*(parameters[parameters.size()-2]*x+parameters.back()); 794 if (step_info) 795 gprintf("Homogeneous solution %gen",makevecteur(sol),step_info,contextptr); 796 bool b=calc_mode(contextptr)==1; 797 if (b) 798 calc_mode(0,contextptr); 799 gen part=_integrate(makesequence(-cst/v.front()*exp(-rac.front()*x,contextptr),x),contextptr)*x+_integrate(makesequence(cst/v.front()*x*exp(-rac.front()*x,contextptr),x),contextptr); 800 if (step_info) 801 gprintf("Particuliar solution %gen",makevecteur(part),step_info,contextptr); 802 if (b) 803 calc_mode(1,contextptr); 804 part=simplify(part*exp(rac.front()*x,contextptr),contextptr); 805 result=sol+part; 806 if (step_info) 807 gprintf("General solution %gen",makevecteur(result),step_info,contextptr); 808 return true; 809 } 810 if (int(rac.size())==n){ 811 gen sol; bool reel=true; 812 for (int j=0;j<n;){ 813 if (j<n-1 && is_zero(ratnormal(rac[j]-conj(rac[j+1],contextptr),contextptr),contextptr)){ 814 gen racr,raci; 815 reim(rac[j],racr,raci,contextptr); 816 if (is_strictly_positive(-raci,contextptr)) 817 raci=-raci; 818 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 819 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 820 sol += exp(racr*x,contextptr)*(parameters[parameters.size()-2]*cos(raci*x,contextptr)+parameters[parameters.size()-1]*sin(raci*x,contextptr)); 821 j+=2; 822 continue; 823 } 824 if (reel && !is_zero(im(rac[j],contextptr))) 825 reel=false; 826 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 827 sol += parameters.back()*exp(rac[j]*x,contextptr); 828 j++; 829 } 830 if (step_info) 831 gprintf("Homogeneous solution %gen",makevecteur(sol),step_info,contextptr); 832 if (derive(cst,x,contextptr)==0 && !is_zero(v.back())){ 833 result=sol-cst/v.back(); 834 return true; 835 } 836 // variation des constantes 837 gen M_=_vandermonde(rac,contextptr),part=0; 838 if (ckmatrix(M_)){ 839 matrice M=*M_._VECTptr; 840 vecteur c(n); 841 c[n-1]=-_trig2exp(cst,contextptr)/v.front(); 842 c=linsolve(mtran(M),c,contextptr); 843 for (unsigned i=0;i<c.size();++i){ 844 bool b=calc_mode(contextptr)==1; 845 if (b) 846 calc_mode(0,contextptr); 847 gen tmp=_lin(c[i]*exp(-rac[i]*x,contextptr),contextptr); 848 tmp = _integrate(makesequence(tmp,x),contextptr); 849 part += _lin(tmp*exp(rac[i]*x,contextptr),contextptr); 850 if (b) 851 calc_mode(1,contextptr); 852 } 853 if (reel && is_zero(im(cst,contextptr)) && lop(part,at_integrate).empty()) 854 part=re(part,contextptr); 855 if (1) // desolve((y''+y=sin(x)) and (y(0)=1) and (y'(0)=2),y) 856 part=recursive_ratnormal(part,contextptr); 857 else 858 part=simplify(part,contextptr); 859 } 860 if (step_info) 861 gprintf("Particuliar solution %gen",makevecteur(part),step_info,contextptr); 862 result=sol+part; 863 return true; 864 } 865 } // end n<=3 866 gen laplace_cst=_laplace(makesequence(-cst,x,t),contextptr); 867 if (!is_undef(laplace_cst)){ 868 vecteur lopei=mergevecteur(lop(laplace_cst,at_Ei),lop(laplace_cst,at_integrate)); 869 if (lopei.empty()){ 870 gen arbitrary,tmp; 871 for (int i=n-1;i>=0;--i){ 872 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 873 tmp=tmp*t+parameters.back(); 874 arbitrary=arbitrary+v[i]*tmp; 875 } 876 arbitrary=(laplace_cst+arbitrary)/symb_horner(v,t); 877 arbitrary=ilaplace(arbitrary,t,x,contextptr); 878 result=arbitrary; 879 return true; 880 } 881 } 882 } 883 if (n==2){ // a(x)*y''+b(x)*y'+c(x)*y+d(x)=0 884 gen & a=v[0]; 885 gen & b=v[1]; 886 gen & c=v[2]; 887 gen & d=cst; 888 #if 0 889 if (is_exactly_zero(c)){ 890 vecteur v1(makevecteur(a,b,d)); 891 if (desolve_linn(x,y,t,1,v1,parameters,result,step_info,contextptr)){ 892 result=_integrate(makesequence(result,x),contextptr); 893 return true; 894 } 895 } 896 #endif 897 gen u=-b/a,V=-c/a,w=-d/a, 898 k=simplify(u*u/4-derive(u,x,contextptr)/2+V,contextptr); 899 // y''=u*y'+V*y+w (with u,V,w functions of x) 900 // Pseudo-code from fhub on HP Museum Forum 901 /* 902 k:=u^2/4-u'/2+V 903 if k==const or k*x^2=const then 904 if k=const 905 then s:=x; t:=e^(int(u,x)/2); 906 else u:=u*x+1; k:=u^2/4+V*x^2; s:=ln(x); t:=x^(u/2); 907 endif; 908 if k=0 then u:=t*s; V:=t; 909 elseif k>0 then u:=t*e^(sqrt(k)*s); V:=t*e^(-sqrt(k)*s); 910 else u:=t*cos(sqrt(-k)*s); V:=t*sin(sqrt(-k)*s); 911 endif; 912 w:=w/(u*V'-V*u'); w:=V*int(u*w,x)-u*int(V*w,x); 913 solution: y=c1*u+c2*V+w 914 endif 915 */ 916 bool cst=is_zero(derive(k,x,contextptr)); 917 bool x2=is_zero(derive(ratnormal(u*x,contextptr),x,contextptr)) && is_zero(derive(ratnormal(v*x*x,contextptr),x,contextptr)); 918 if (cst || x2){ 919 gen s,t; 920 if (cst){ 921 s=x; 922 t=simplify(exp(integrate_without_lnabs(u,x,contextptr)/2,contextptr),contextptr); 923 } 924 else { 925 u=u*x+1; 926 u=simplify(u,contextptr); 927 k=simplify(u*u/4+V*x*x,contextptr); 928 s=ln(x,contextptr); t=pow(x,u/2,contextptr); 929 } 930 if (is_zero(k)){ 931 u=t*s; V=t; 932 } 933 else { 934 if (is_strictly_positive(-k,contextptr)){ 935 gen tmp=sqrt(-k,contextptr)*s; 936 u=t*cos(tmp,contextptr); 937 V=t*sin(tmp,contextptr); 938 } 939 else { 940 if (s.is_symb_of_sommet(at_ln)){ 941 gen tmp=pow(s._SYMBptr->feuille,sqrt(k,contextptr),contextptr); 942 u=t*tmp; 943 V=t/tmp; 944 } 945 else { 946 gen tmp=sqrt(k,contextptr)*s; 947 u=t*exp(tmp,contextptr); 948 V=t*exp(-tmp,contextptr); 949 } 950 } 951 } 952 w=simplify(w/(u*derive(V,x,contextptr)-V*derive(u,x,contextptr)),contextptr); 953 w=V*integrate_without_lnabs(u*w,x,contextptr)- 954 u*integrate_without_lnabs(V*w,x,contextptr); 955 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 956 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 957 result=w+parameters[parameters.size()-2]*u+parameters[parameters.size()-1]*V; 958 return true; 959 } 960 // IMPROVE: if a, b, c are polynomials, search for a polynomial solution 961 // of the homogeneous equation, if found we can solve the diffeq 962 gen aa(a),bb(b),cc(c); 963 if (lvarxwithinv(makevecteur(a,b,c),x,contextptr)==vecteur(1,x)){ 964 vecteur l=vecteur(1,x); 965 gen a0(a),b0(b); 966 a=_coeff(makesequence(a,x),contextptr); 967 b=_coeff(makesequence(b,x),contextptr); 968 c=_coeff(makesequence(c,x),contextptr); 969 if (a.type==_VECT && b.type==_VECT && c.type==_VECT){ 970 int A=int(a._VECTptr->size())-1,B=int(b._VECTptr->size())-1,C=int(c._VECTptr->size())-1,N=-1; 971 if (C==B-1){ 972 gen n=-c._VECTptr->front()/b._VECTptr->front(); 973 if (n.type==_INT_ && n.val>N){ 974 if (A-2<C || n==1) 975 N=n.val; 976 } 977 if (A-2==C){ 978 // a*n*(n-1)+b*n+c=a*n^2+(b-1)*n+c=0 979 gen aa=a._VECTptr->front(),bb=b._VECTptr->front()-1,cc=c._VECTptr->front(); 980 gen delta=(sqrt(bb*bb-4*aa*cc,contextptr)+bb)/2; 981 if (delta.type==_INT_ && delta.val>N) 982 N=delta.val; 983 } 984 } 985 if (A-2==B-1 && C<B-1){ 986 gen n=-b._VECTptr->front()/a._VECTptr->front()+1; 987 if (n.type==_INT_ && n.val>N) 988 N=n.val; 989 } 990 if (C==A-2 && B-1<C){ 991 gen delta=(1+sqrt(1+4*c._VECTptr->front()/a._VECTptr->front(),contextptr))/2; 992 if (delta.type==_INT_ && delta.val>N) 993 N=delta.val; 994 } 995 if (N>=0){ 996 int nrows=giacmax(giacmax(B,C+1),N==1?0:A)+N; 997 // search a solution sum(y_k*x*k,k,0,N) 998 matrice m(nrows); 999 for (int i=0;i<nrows;++i) 1000 m[i]=vecteur(N+1); 1001 // a*y'' 1002 for (int i=0;i<a._VECTptr->size();++i){ 1003 int j=int(a._VECTptr->size())-i-1; 1004 for (int k=2;k<=N;++k){ 1005 (*m[j+k-2]._VECTptr)[k] += k*(k-1)*a[i]; 1006 } 1007 } 1008 // b*y' 1009 for (int i=0;i<b._VECTptr->size();++i){ 1010 int j=int(b._VECTptr->size())-i-1; 1011 for (int k=1;k<=N;++k){ 1012 (*m[j+k-1]._VECTptr)[k] += k*b[i]; 1013 } 1014 } 1015 // c*y 1016 for (int i=0;i<c._VECTptr->size();++i){ 1017 int j=int(c._VECTptr->size())-i-1; 1018 for (int k=0;k<=N;++k){ 1019 (*m[j+k]._VECTptr)[k] += c[i]; 1020 } 1021 } 1022 m=mker(m,contextptr); 1023 if (!m.empty()){ 1024 gen sol=m.front(); 1025 if (sol.type==_VECT){ 1026 vecteur v=*sol._VECTptr; 1027 reverse(v.begin(),v.end()); 1028 sol=symb_horner(-v,x); 1029 *logptr(contextptr) << "Polynomial solution found " << sol << '\n'; 1030 // now solve equation a*y''+b*y'+c*y+d=0 with y=sol*z 1031 // a*sol*z''+(2*a*sol'+b*sol)*z'=d 1032 gen res=desolve_lin1(a0*sol,2*a0*derive(sol,x,contextptr)+b0*sol,d,x,parameters,step_info,contextptr); 1033 res=_integrate(makesequence(res,x),contextptr); 1034 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 1035 res += parameters.back(); 1036 res=res*sol; 1037 result=res; 1038 return true; 1039 } 1040 } 1041 } 1042 } 1043 } // end polynomial a,b,c 1044 #ifndef USE_GMP_REPLACEMENTS 1045 a=aa; b=bb; c=cc; 1046 if (d==0 && lvarx(makevecteur(a,b,c),x,contextptr)==vecteur(1,x)){ 1047 // if a,b,c are rationals and d==0, Kovacic 1048 gen k=_kovacicsols(makesequence(makevecteur(a,b,c),x),contextptr); 1049 if (k.type==_VECT && !k._VECTptr->empty()){ 1050 if (k._VECTptr->size()==2){ 1051 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 1052 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 1053 result=parameters[parameters.size()-2]*k._VECTptr->front()+parameters[parameters.size()-1]*k._VECTptr->back(); 1054 return true; 1055 } 1056 if (k._VECTptr->size()==1){ 1057 // we have one solution Y, find an independent one as z*Y 1058 gen Y=k._VECTptr->front(); 1059 // a*(zY)''+b*(zY)'+c*(zY)=0 1060 // a*(z''Y+2z'Y')+b*(z'Y)=0 1061 // z''*(aY)+z'*(2aY'+bY)=0 1062 result=desolve_lin1(a*Y,2*a*derive(Y,x,contextptr)+b*Y,0,x,parameters,step_info,contextptr); 1063 result=_integrate(makesequence(result,x),contextptr); 1064 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 1065 result += parameters.back(); 1066 result = result*Y; 1067 return true; 1068 } 1069 } 1070 } 1071 #endif 1072 } // end 2nd order eqdiff 1073 return false; 1074 } 1075 desolve_f(const gen & f_orig,const gen & x_orig,const gen & y_orig,int & ordre,vecteur & parameters,gen & fres,int step_info,bool & num,GIAC_CONTEXT)1076 gen desolve_f(const gen & f_orig,const gen & x_orig,const gen & y_orig,int & ordre,vecteur & parameters,gen & fres,int step_info,bool & num,GIAC_CONTEXT){ 1077 num=false; 1078 // if x_orig.type==_VECT || y_orig.type==_VECT, they should be evaled 1079 if (x_orig.type!=_VECT && eval(x_orig,1,contextptr)!=x_orig) 1080 return gensizeerr("Independant variable assigned. Run purge("+x_orig.print(contextptr)+")\n"); 1081 if (y_orig.type!=_VECT && eval(y_orig,1,contextptr)!=y_orig) 1082 return gensizeerr("Dependant variable assigned. Run purge("+y_orig.print(contextptr)+")\n"); 1083 gen x(x_orig); 1084 if ( (x_orig.type==_VECT) && (x_orig._VECTptr->size()==1) ) 1085 x=x_orig._VECTptr->front(); 1086 if (x.type!=_IDNT){ 1087 gen vx,vy; 1088 ggb_varxy(f_orig,vx,vy,contextptr); 1089 if (x_orig.type==_VECT) 1090 return desolve_with_conditions(makevecteur(f_orig,x_orig,y_orig),vx,vy,fres,step_info,contextptr); 1091 else 1092 return desolve_with_conditions(makevecteur(f_orig,makevecteur(x_orig,y_orig)),vx,vy,fres,step_info,contextptr); 1093 } 1094 if (y_orig.type==_VECT) // FIXME: differential system 1095 return gensizeerr(contextptr); 1096 gen f=remove_and(f_orig,at_and); 1097 if (f.type==_VECT){ 1098 vecteur fv=*f._VECTptr; 1099 return desolve_with_conditions(fv,x,y_orig,fres,step_info,contextptr); 1100 } 1101 gen y(y_orig),yof(y_orig),partic(undef); 1102 if (y_orig.is_symb_of_sommet(at_equal)){ 1103 // particular solution provided 1104 y=y_orig._SYMBptr->feuille[0]; 1105 partic=eval(y_orig._SYMBptr->feuille[1],1,contextptr); 1106 } 1107 if (y.type==_IDNT){ 1108 yof=symb_of(y,gen(vecteur(1,x),_SEQ__VECT)); 1109 f=quotesubst(f,yof,y,contextptr); 1110 f=quotesubst(f,y,yof,contextptr); 1111 } 1112 else 1113 y=function_of(y_orig,x); 1114 if (is_undef(y)) 1115 return y; 1116 gen save_vx=vx_var; 1117 vx_var=x; 1118 int save=calc_mode(contextptr); 1119 calc_mode(0,contextptr); 1120 #ifdef GIAC_HAS_STO_38 1121 // HP Prime: if there is a M0-M9 identifier, this will not work 1122 vecteur fid(lidnt(f)); 1123 for (unsigned i=0;i<fid.size();++i){ 1124 gen fidi=fid[i]; 1125 if (fidi.type==_IDNT){ 1126 const char * ch=fidi._IDNTptr->id_name; 1127 if (strlen(ch)==2 && ch[0]=='M' && ch[1]>='0' && ch[1]<='9') 1128 return gensizeerr("Home matrix variable "+ string(ch)+" not allowed. Store your matrix in a CAS variable first"); 1129 } 1130 } 1131 #endif 1132 f=remove_equal(eval(f,eval_level(contextptr),contextptr)); 1133 num=has_num_coeff(f); 1134 if (num) 1135 f=exact(f,contextptr); 1136 if (ckmatrix(f)){ 1137 vecteur v = *f._VECTptr; 1138 for (int i=0;i<v.size();++i){ 1139 v[i].subtype=0; 1140 } 1141 f=v; 1142 } 1143 calc_mode(save,contextptr); 1144 fres=f=quotesubst(f,yof,y,contextptr); 1145 vx_var=save_vx; 1146 // Here f= f(derive(y,x),y) for a 1st order equation 1147 int n=diffeq_order(f,y); 1148 if (n==0) 1149 return solve(f,y,0,contextptr); 1150 if (n<=0) 1151 return gensizeerr(contextptr); 1152 vecteur v; 1153 gen t=gen_t(makevecteur(x,y,f),contextptr); 1154 if (is_linear_diffeq(f,x,y,n,v,step_info,contextptr)){ 1155 gen result; 1156 if (n>1 && !is_undef(partic)){ 1157 // reduce order by one 1158 vecteur s(n,partic); 1159 for (int i=1;i<n;++i){ 1160 s[i]=derive(s[i-1],x,contextptr); 1161 } 1162 vecteur w(n+1); 1163 w[n]=v[n+1]; // cst coeff 1164 for (int l=0;l<n;++l){ 1165 gen tmp=0; 1166 for (int j=0;j<=l;++j){ 1167 tmp += v[j]*comb(n-j,l-j)*s[l-j]; 1168 } 1169 w[l]=tmp; 1170 } 1171 if (desolve_linn(x,y,t,n-1,w,parameters,result,step_info,contextptr)){ 1172 result=integrate_without_lnabs(result,x,contextptr); 1173 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 1174 result = partic*(result+parameters.back()); 1175 return result; 1176 } 1177 } 1178 if (desolve_linn(x,y,t,n,v,parameters,result,step_info,contextptr)) 1179 return result; 1180 } 1181 vecteur substin(n); 1182 vecteur substout(n); 1183 for (int i=0;i<n;++i){ 1184 substin[i]=symb_derive(y,x,i+1); 1185 substout[i]=identificateur(" y"+print_INT_(i)); 1186 } 1187 gen ff=quotesubst(f,substin,substout,contextptr); 1188 if (is_zero(derive(ff,y,contextptr))){ // y incomplete 1189 if (step_info) 1190 gprintf("y-incomplete",vecteur(0),step_info,contextptr); 1191 for (int i=0;i<n;++i){ 1192 substout[i]=symb_derive(y,x,i); 1193 } 1194 f=quotesubst(f,substin,substout,contextptr); 1195 int tmp; 1196 gen sol=desolve(f,x,y,tmp,parameters,contextptr); 1197 if (is_undef(sol)) return sol; 1198 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 1199 gen p(parameters.back()); 1200 if (sol.type==_VECT) 1201 p=vecteur(sol._VECTptr->size(),p); 1202 sol=integrate_without_lnabs(sol,x,contextptr)+p; 1203 return sol; 1204 } 1205 if (n==1) { // 1st order 1206 vecteur sol; 1207 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 1208 f=quotesubst(f,symb_derive(y,x),t,contextptr); 1209 // f is an expression of x,y,t where t stands for y' 1210 gen fa,fb,fc,fd,faa,fab; 1211 // Test for Lagrange/Clairault-like eqdiff, 1212 if (x.type==_IDNT && y.type==_IDNT && is_linear_wrt(f,y,fc,fd,contextptr) && is_linear_wrt(fd,x,fa,fb,contextptr)){ 1213 // Clairault: fa must be cst*t and fc must be cst (must simplify fa and fc) 1214 // f=y*fc+(fa*x+fb) 1215 fd=gcd(fc,fa); 1216 fa=normal(fa/fd,contextptr); fb=normal(fb/fd,contextptr); fc=normal(fc/fd,contextptr); 1217 if (is_linear_wrt(fa,t,faa,fab,contextptr) && is_zero(fab) && derive(faa,makevecteur(x,y,t),contextptr)==vecteur(3,0) && derive(fc,makevecteur(x,y,t),contextptr)==vecteur(3,0) && derive(fb,makevecteur(x,y),contextptr)==vecteur(2,0)){ 1218 // 0=f=fc*y+fd = fc*y+fa*x+fb = fc*y+faa*x*y'+fb 1219 // -> y=-faa/fc*x*y' -fb/fc 1220 if (is_one(ratnormal(-faa/fc,contextptr))){ 1221 if (step_info) 1222 gprintf("Order 1 Clairault differential equation",vecteur(0),step_info,contextptr); 1223 // y=x*y'-fb/fc 1224 gen fm=ratnormal(-fb/fc,contextptr); 1225 gen fmp=derive(fm,t,contextptr); 1226 sol.push_back(parameters.back()*x+subst(fm,t,parameters.back(),false,contextptr)); 1227 sol.push_back(makevecteur(-fmp,-t*fmp+fm)); 1228 return sol; 1229 } 1230 } 1231 // Lagrange-> fa/fb/fc dependent de t uniquement, if fb==0 -> separate var or homogeneous 1232 if (is_zero(derive(makevecteur(fa,fb,fc),x,contextptr)) && !is_zero(fb)){ 1233 if (step_info) 1234 gprintf("Order 1 Lagrange differential equation",vecteur(0),step_info,contextptr); 1235 // y+fa/fc*x+fb/fc=0 1236 fa=fa/fc; fb=fb/fc; 1237 // y+fa*x+fb=0 1238 // t=dy/dx, dy/dt=t*dx/dt => t*dx/dt+fa'*x+fb'+fa*dx/dt 1239 // linear equation 1st order (fa+t)*dx/dt+fa'*x+fb'=0 1240 gen res=desolve_lin1(fa+t,derive(fa,t,contextptr),derive(fb,t,contextptr),t,parameters,step_info,contextptr); 1241 vecteur sing(solve(t+fa,t,3,contextptr)); 1242 for (int i=0;i<int(sing.size());++i){ 1243 sing[i]=subst(-fa*x-fb,t,sing[i],false,contextptr); 1244 } 1245 // should deparametrize like for homogeneous if possible 1246 #ifdef NO_STDEXCEPT 1247 vecteur newsol=solve(res-x,*t._IDNTptr,3,contextptr); 1248 if (is_undef(newsol)){ 1249 newsol.clear(); 1250 *logptr(contextptr) << "Unable to solve implicit equation "<< res-x << "=0 in " << t << '\n'; 1251 } 1252 #else 1253 vecteur newsol; 1254 try { 1255 newsol=solve(res-x,*t._IDNTptr,3,contextptr); 1256 } catch(std::runtime_error & err){ 1257 last_evaled_argptr(contextptr)=NULL; 1258 newsol.clear(); 1259 *logptr(contextptr) << "Unable to solve implicit equation "<< res-x << "=0 in " << t << '\n'; 1260 } 1261 #endif 1262 if (newsol.empty()) 1263 sing.push_back(makevecteur(res,-fa*res-fb)); 1264 else { 1265 for (int i=0;i<int(newsol.size());++i){ 1266 sing.push_back(subst(-fa*x-fb,t,newsol[i],false,contextptr)); 1267 } 1268 } 1269 return sing; 1270 } 1271 } // end Lagrange-Clairault 1272 vecteur v(solve(f,t,3,contextptr)); // now solve y'=v[i](y) 1273 const_iterateur it=v.begin(),itend=v.end(); 1274 for (;it!=itend;++it){ 1275 // Separable variables? 1276 f=factors(*it,x,contextptr); // Factor then split factors 1277 gen xfact(plus_one),yfact(plus_one); 1278 if (separate_variables(f,x,y,xfact,yfact,step_info,contextptr)){ // y'/yfact=xfact 1279 gen pr=integrate_without_lnabs(inv(yfact,contextptr),y,contextptr); 1280 #if 1 1281 vecteur prv=lop(lvarx(pr,y),at_ln); 1282 gen pra,prb; 1283 if (!prv.empty() && prv[0].is_symb_of_sommet(at_ln) && is_linear_wrt(pr,prv[0],pra,prb,contextptr)){ 1284 pr=_lncollect(pra*(symbolic(at_ln,parameters.back()*prv[0]._SYMBptr->feuille))+prb,contextptr); 1285 } 1286 else 1287 pr=parameters.back()+pr; 1288 #else 1289 if (has_op(pr,*at_ln)) 1290 pr=_lncollect(pr,contextptr); // hack to solve y'=y*(1-y) 1291 if (pr.is_symb_of_sommet(at_ln)) 1292 pr=symbolic(at_ln,parameters.back()*pr._SYMBptr->feuille); 1293 else 1294 pr=parameters.back()+pr; 1295 #endif 1296 gen implicitsol=pr-integrate_without_lnabs(xfact,x,contextptr); 1297 #ifdef NO_STDEXCEPT 1298 vecteur newsol=solve(implicitsol,*y._IDNTptr,3,contextptr); 1299 if (is_undef(newsol)){ 1300 newsol.clear(); 1301 *logptr(contextptr) << "Unable to solve implicit equation "<< implicitsol << "=0 in " << y << '\n'; 1302 } 1303 #else 1304 vecteur newsol; 1305 int cm=calc_mode(contextptr); 1306 calc_mode(0,contextptr); 1307 try { 1308 newsol=solve(implicitsol,*y._IDNTptr,3,contextptr); 1309 } catch(std::runtime_error & err){ 1310 last_evaled_argptr(contextptr)=NULL; 1311 newsol.clear(); 1312 *logptr(contextptr) << "Unable to solve implicit equation "<< implicitsol << "=0 in " << y << '\n'; 1313 } 1314 calc_mode(cm,contextptr); 1315 #endif 1316 sol=mergevecteur(sol,newsol); 1317 continue; 1318 } // end separate variables 1319 if (is_zero(derive(*it,x,contextptr))){ // x incomplete 1320 if (step_info) 1321 gprintf("Order 1 x-incomplete differential equation",vecteur(0),step_info,contextptr); 1322 if (debug_infolevel) 1323 *logptr(contextptr) << gettext("Incomplete") << '\n'; 1324 gen pr=integrate_without_lnabs(inv(*it,contextptr),y,contextptr)+parameters.back(); 1325 sol=mergevecteur(sol,solve(pr-x,*y._IDNTptr,3,contextptr)); 1326 continue; 1327 } 1328 // check for a linear substitution -> like an x incomplete 1329 fa=derive(*it,x,contextptr); fb=derive(*it,y,contextptr); 1330 fc=simplify(fa/fb,contextptr); 1331 if (is_zero(derive(fc,x,contextptr)) && is_zero(derive(fc,y,contextptr))){ 1332 gen eff=subst(*it,y,y-fc*x,false,contextptr); // does not depend on x 1333 gen pr=integrate_without_lnabs(inv(eff+fc,contextptr),y,contextptr)+parameters.back(); 1334 pr=subst(pr,y,y+fc*x,false,contextptr); 1335 vecteur l1=lop(lvarx(pr,y),at_floor); 1336 if (!l1.empty()){ 1337 vecteur l2(l1.size()); 1338 pr=subst(pr,l1,l2,false,contextptr); 1339 } 1340 vecteur sol1=solve(pr-x,*y._IDNTptr,3,contextptr); 1341 sol=mergevecteur(sol,sol1); 1342 continue; 1343 } 1344 // homogeneous? 1345 gen tplus(t); 1346 gen tmpsto=sto(doubleassume_and(vecteur(2,0),0,1,false,contextptr),tplus,contextptr); 1347 if (is_undef(tmpsto)) 1348 return tmpsto; 1349 f=quotesubst(*it,makevecteur(x,y),makevecteur(tplus*x,tplus*y),contextptr); 1350 f=recursive_normal(f-*it,contextptr); 1351 purgenoassume(tplus,contextptr); 1352 if (is_zero(f)){ 1353 if (step_info) 1354 gprintf("Order 1 Homogeneous differential equation",vecteur(0),step_info,contextptr); 1355 if (debug_infolevel) 1356 *logptr(contextptr) << gettext("Homogeneous differential equation") << '\n'; 1357 tmpsto=sto(doubleassume_and(vecteur(2,0),0,1,false,contextptr),x,contextptr); 1358 if (is_undef(tmpsto)) 1359 return tmpsto; 1360 f=recursive_normal(quotesubst(*it,y,tplus*x,contextptr)-tplus,contextptr); 1361 purgenoassume(x,contextptr); 1362 // y=tx -> t'x=f 1363 // Singular solutions f(t)=0 1364 vecteur singuliere(multvecteur(x,solve(f,t,complex_mode(contextptr) + 2,contextptr))); 1365 sol=mergevecteur(sol,singuliere); 1366 // Non singular: t'/f(t)=1/x 1367 gen pr=parameters.back()*_simplify(exp(integrate_without_lnabs(inv(f,contextptr),t,contextptr),contextptr),contextptr); 1368 // Try to find t in x=pr 1369 vecteur v=protect_solve(x-pr,*t._IDNTptr,1,contextptr); 1370 if (!v.empty() && !is_undef(v)){ 1371 *logptr(contextptr) << "solve(" << pr << "=" << x << "," << t << ") returned " << v << ".\nIf solutions were missed consider paramplot(" << makevecteur(pr,t*pr) << "," << t << ")" << '\n'; 1372 for (unsigned j=0;j<v.size();++j){ 1373 sol.push_back(x*v[j]); 1374 } 1375 } 1376 else 1377 sol.push_back(gen(makevecteur(pr,t*pr),_CURVE__VECT)); 1378 continue; 1379 } 1380 // exact? y'=*it=f(x,y) -> N dy + M dx=0 where -M/N=y' 1381 gen M,N; 1382 f=_fxnd(*it,contextptr); 1383 M=-f[0]; 1384 N=f[1]; 1385 // find an integrating factor P such that d_x(P*N)=d_y(P*M) 1386 // If P depends on x then N*d_x(P)+Pd_x(N)=Pd_y(M) -> 1387 // d_x(P)/P=(d_y(M)-d_x(N))/N should depend on x only 1388 // If P depends on y then P d_x(N)=Pd_y(M)+Md_y(P) 1389 // d_y(P)/P=(d_x(N)-d_y(M))/M 1390 // Then solve P*Ndy+P*Mdx=dF 1391 f=normal((derive(M,y,contextptr)-derive(N,x,contextptr))/N,contextptr); 1392 if (is_zero(derive(f,y,contextptr))){ 1393 gen P=simplify(exp(integrate_without_lnabs(f,x,contextptr),contextptr),contextptr); 1394 // D_y(F)=P*N 1395 gen F=P*integrate_without_lnabs(N,y,contextptr); 1396 if (step_info) 1397 gprintf("Order 1 Integrating factor %gen",makevecteur(P),step_info,contextptr); 1398 // D_x(F)=P*M 1399 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 1400 F=F+integrate_without_lnabs(normal(P*M-derive(F,x,contextptr),contextptr),x,contextptr)+parameters.back(); 1401 sol=mergevecteur(sol,solve(F,*y._IDNTptr,3,contextptr)); 1402 continue; 1403 } 1404 f=normal((derive(N,x,contextptr)-derive(M,y,contextptr))/M,contextptr); 1405 if (is_zero(derive(f,x,contextptr))){ 1406 gen P=simplify(exp(integrate_without_lnabs(f,y,contextptr),contextptr),contextptr); 1407 gen F=P*integrate_without_lnabs(M,x,contextptr); 1408 // D_y(F)=P*N 1409 if (step_info) 1410 gprintf("Order 1 Integrating factor %gen",makevecteur(P),step_info,contextptr); 1411 F=F+integrate_without_lnabs(normal(P*N-derive(F,y,contextptr),contextptr),y,contextptr)+diffeq_constante(int(parameters.size()),contextptr); 1412 sol=mergevecteur(sol,solve(F,*y._IDNTptr,3,contextptr)); 1413 continue; 1414 } 1415 // Bernoulli? 1416 // y'=a(x)*y+b(x)*y^k 1417 // Let z=y^(1-k) 1418 // z'=(1-k)*y^(-k)*y'=(1-k)*[a(x)*z+b(x)] 1419 // Solve for z then for y 1420 f=subst(*it,y,2*y,false,contextptr); 1421 f=factors(f-2*(*it),vx_var,contextptr); // should be (2^k-2)*b(x)*y^k 1422 xfact=plus_one; 1423 yfact=plus_one; 1424 if (separate_variables(f,x,y,xfact,yfact,step_info,contextptr)){ 1425 // xfact should be (2^k-2)*b(x) and yfact=y^k 1426 if ( (yfact.type==_SYMB) && (yfact._SYMBptr->sommet==at_pow) && 1427 (yfact._SYMBptr->feuille._VECTptr->front()==y) ){ 1428 if (step_info) 1429 gprintf("Order 1 Bernoulli differential equation",vecteur(0),step_info,contextptr); 1430 gen k=yfact._SYMBptr->feuille._VECTptr->back(); 1431 gen B=normal(xfact/(pow(plus_two,k,contextptr)-plus_two),contextptr); 1432 gen A=normal((*it-B*pow(y,k,contextptr))/y,contextptr); 1433 gen b=(k-1)*A; 1434 gen c=(k-1)*B; 1435 gen i=simplify(integrate_without_lnabs(b,x,contextptr),contextptr); 1436 gen C=integrate_without_lnabs(-c*exp(i,contextptr),x,contextptr); 1437 f= (C+parameters.back())*exp(-i,contextptr); 1438 gen sol1=pow(f,inv(1-k,contextptr),contextptr); 1439 sol.push_back(sol1); 1440 // FIXME: we should add other roots of unity in complex mode 1441 if (k.type==_INT_ && k.val %2) 1442 sol.push_back(-sol1); 1443 } 1444 } 1445 // Ricatti f=*it quadratic in y 1446 gen P,Q,R; 1447 if (is_quadratic_wrt(*it,y,R,Q,P,contextptr)){ 1448 if (step_info) 1449 gprintf("Order 1 Riccati differential equation",vecteur(0),step_info,contextptr); 1450 gen result; 1451 // y'=P+Q*y+R*y^2=q0+q1*y+q2*y^2 1452 if (!is_undef(partic)){ 1453 // z'+(q1+2*q2*partic)*z+q2=0 1454 result=desolve_lin1(1,Q+2*R*partic,R,x,parameters,step_info,contextptr); 1455 return makevecteur(partic,partic+inv(result,contextptr)); 1456 } 1457 // let y=-1/(R*F)*dF/dx, then F''-(1/R*R'+Q)*F'+R*P*F=0 1458 vecteur v(makevecteur(1,-normal(Q+derive(R,x,contextptr)/R,contextptr),normal(R*P,contextptr),0)); 1459 if (desolve_linn(x,y,t,2,v,parameters,result,step_info,contextptr)){ 1460 result=lnexpand(ln(result,contextptr),contextptr); 1461 result=-derive(result,x,contextptr)/R; 1462 result=ratnormal(result,contextptr); 1463 gen lastp=parameters.back(); 1464 parameters.pop_back(); 1465 gen partic=subst(result,lastp,0,false,contextptr); 1466 partic=ratnormal(partic,contextptr); 1467 result=subst(result,lastp,1,false,contextptr); 1468 result=ratnormal(result,contextptr); 1469 //result=-derive(result,x,contextptr)/(R*result); 1470 return makevecteur(partic,result); 1471 } 1472 } 1473 } // end for (;it!=itend;) 1474 return sol; 1475 } // end if n==1 1476 if (n==2){ 1477 // y''=f(y,y'), set u=y' -> u'=f(y,u)/u 1478 gen der1=substout[0],der2=substout[1]; 1479 gen soly2=_cSolve(makesequence(symb_equal(ff,0),der2),contextptr); 1480 vecteur paramsave=parameters; 1481 if (soly2.type==_VECT && !is_undef(soly2)){ 1482 vecteur sol; 1483 const vecteur & soly2v = *soly2._VECTptr; 1484 for (unsigned i=0;i<soly2v.size();++i){ 1485 gen soly2c=soly2v[i]; 1486 gen a,b,c; 1487 if (is_quadratic_wrt(soly2c,der1,a,b,c,contextptr) 1488 && is_zero(c) && is_zero(derive(a,x,contextptr)) 1489 && is_zero(derive(b,y,contextptr)) ){ 1490 parameters=paramsave; 1491 parameters.push_back(diffeq_constante(int(parameters.size()),contextptr)); 1492 gen usolj=parameters.back()*exp(integrate_without_lnabs(b,x,contextptr),contextptr)*exp(integrate_without_lnabs(a,y,contextptr),contextptr); 1493 gen ysol=desolve(symb_equal(symbolic(at_derive,makesequence(y,x)),usolj),x,y,ordre,parameters,contextptr); 1494 if (is_undef(ysol)) 1495 return unable_to_solve_diffeq(); 1496 sol=mergevecteur(sol,gen2vecteur(ysol)); 1497 continue; 1498 } 1499 if (is_zero(derive(soly2c,x,contextptr))){ // x-incomplete 1500 if (step_info) 1501 gprintf("Order 2 x-incomplete differential equation",vecteur(0),step_info,contextptr); 1502 // desolve(u'=soly2c/der1,y,u) 1503 parameters=paramsave; 1504 gen usol=desolve(symb_equal(symbolic(at_derive,makesequence(der1,y)),soly2c/der1),y,der1,ordre,parameters,contextptr); 1505 if (is_undef(usol)) 1506 return unable_to_solve_diffeq(); 1507 if (usol.type!=_VECT) 1508 usol=vecteur(1,usol); 1509 vecteur paramsavein=parameters; 1510 for (unsigned j=0;j<usol._VECTptr->size();++j){ 1511 parameters=paramsavein; 1512 gen usolj=(*usol._VECTptr)[j]; 1513 gen ysol=desolve(symb_equal(symbolic(at_derive,makesequence(y,x)),usolj),x,y,ordre,parameters,contextptr); 1514 if (is_undef(ysol)) 1515 return unable_to_solve_diffeq(); 1516 sol=mergevecteur(sol,gen2vecteur(ysol)); 1517 } 1518 continue; 1519 } // end x-incomplete 1520 gen res(string2gen(gettext("Unable to solve differential equation"),false)); 1521 res.subtype=-1; 1522 sol.push_back(res); 1523 } 1524 ordre=2; 1525 return sol; 1526 } 1527 } 1528 return unable_to_solve_diffeq(); 1529 } ggbputinlist(const gen & g,GIAC_CONTEXT)1530 gen ggbputinlist(const gen & g,GIAC_CONTEXT){ 1531 if (g.type==_VECT || calc_mode(contextptr)!=1) 1532 return g; 1533 return makevecteur(g); 1534 } point2vecteur(const gen & g_,GIAC_CONTEXT)1535 static gen point2vecteur(const gen & g_,GIAC_CONTEXT){ 1536 if (!g_.is_symb_of_sommet(at_point)) 1537 return g_; 1538 gen g=g_._SYMBptr->feuille; 1539 gen x,y; 1540 if (g.type==_VECT){ 1541 if (g._VECTptr->size()!=2) 1542 return gensizeerr(contextptr); 1543 x=g._VECTptr->front(); 1544 y=g._VECTptr->back(); 1545 } 1546 else 1547 reim(g,x,y,contextptr); 1548 g=makevecteur(x,y); 1549 return g; 1550 } 1551 // "unary" version _desolve(const gen & args,GIAC_CONTEXT)1552 gen _desolve(const gen & args,GIAC_CONTEXT){ 1553 if ( args.type==_STRNG && args.subtype==-1) return args; 1554 //if (has_num_coeff(args)) return evalf(_desolve(exact(args,contextptr),contextptr),1,contextptr); 1555 int ordre; 1556 vecteur parameters; 1557 if (args.type!=_VECT || args.subtype!=_SEQ__VECT || (!args._VECTptr->empty() && is_equal(args._VECTptr->back()) && args._VECTptr->back()._SYMBptr->feuille[0].type!=_IDNT)){ 1558 // guess x and y 1559 vecteur lv(lop(args,at_of)); 1560 vecteur f; 1561 if (lv.size()>=1 && lv[0]._SYMBptr->feuille.type==_VECT && (f=*lv[0]._SYMBptr->feuille._VECTptr).size()==2){ 1562 if (f[1].type==_IDNT || f[1].is_symb_of_sommet(at_at)){ 1563 return desolve(args,f[1],f[0],ordre,parameters,contextptr); 1564 } 1565 } 1566 gen vx,vy; 1567 lv=lidnt(evalf(args,1,contextptr)); 1568 if (lv.size()==2){ 1569 vx=lv[0]; 1570 vy=lv[1]; 1571 lv=lvar(apply(args,equal2diff)); 1572 lv=lop(lv,at_derive); 1573 lv=lidnt(lv); 1574 if (lv.size()==1 && vx==lv.front()) 1575 swapgen(vx,vy); 1576 return _desolve(makesequence(args,vx,vy),contextptr); 1577 } 1578 ggb_varxy(args,vx,vy,contextptr); 1579 return _desolve(makesequence(args,vx,vy),contextptr); 1580 } 1581 vecteur v=*args._VECTptr; 1582 int s=int(v.size()); 1583 for (int i=0;i<s;++i){ 1584 v[i]=apply(v[i],point2vecteur,contextptr); 1585 } 1586 if (s==3 && v[1].type==_VECT && v[2].type==_VECT) 1587 swapgen(v[1],v[2]); 1588 if (s==2 && v[1].type==_VECT && v[1]._VECTptr->size()==2){ 1589 gen a=eval(v[1]._VECTptr->front(),1,contextptr); 1590 gen b=eval(v[1]._VECTptr->back(),1,contextptr); 1591 v[1]=a; 1592 v.insert(v.begin()+2,b); 1593 ++s; 1594 } 1595 if (s==2){ 1596 if ( (v[1].type==_SYMB && v[1]._SYMBptr->sommet==at_of && v[1]._SYMBptr->feuille.type==_VECT &&v [1]._SYMBptr->feuille._VECTptr->size()==2 ) ) 1597 return desolve(v[0],(*v[1]._SYMBptr->feuille._VECTptr)[1],(*v[1]._SYMBptr->feuille._VECTptr)[0],ordre,parameters,contextptr); 1598 return ggbputinlist(desolve( v[0],vx_var,v[1],ordre,parameters,contextptr),contextptr); 1599 } 1600 gen f; 1601 if (s==4) 1602 return ggbputinlist(desolve_with_conditions(makevecteur(v[0],v[3]),v[1],v[2],f,contextptr),contextptr); 1603 if (s==5) 1604 return ggbputinlist(desolve_with_conditions(makevecteur(v[0],v[3],v[4]),v[1],v[2],f,contextptr),contextptr); 1605 if (s!=3) 1606 return gensizeerr(contextptr); 1607 return ggbputinlist(desolve( v[0],v[1],v[2],ordre,parameters,contextptr),contextptr); 1608 } 1609 static const char _desolve_s []="desolve"; 1610 static define_unary_function_eval_quoted (__desolve,&_desolve,_desolve_s); 1611 define_unary_function_ptr5( at_desolve ,alias_at_desolve,&__desolve,1,true); 1612 1613 static const char _dsolve_s []="dsolve"; 1614 static define_unary_function_eval_quoted (__dsolve,&_desolve,_dsolve_s); 1615 define_unary_function_ptr5( at_dsolve ,alias_at_dsolve,&__dsolve,_QUOTE_ARGUMENTS,true); 1616 ztrans(const gen & f,const gen & x,const gen & s,GIAC_CONTEXT)1617 gen ztrans(const gen & f,const gen & x,const gen & s,GIAC_CONTEXT){ 1618 if (x.type!=_IDNT) 1619 return gensizeerr(contextptr); 1620 gen t(s); 1621 if (s==x){ 1622 #ifdef GIAC_HAS_STO_38 1623 t=identificateur("z38_"); 1624 #else 1625 t=identificateur(" tztrans"); 1626 #endif 1627 } 1628 if (!assume_t_in_ab(t,plus_inf,plus_inf,true,true,contextptr)) 1629 return gensizeerr(contextptr); 1630 gen tmp=expand(f*pow(t,-x,contextptr),contextptr); 1631 gen res=_sum(gen(makevecteur(tmp,x,0,plus_inf),_SEQ__VECT),contextptr); 1632 purgenoassume(t,contextptr); 1633 if (s==x) 1634 res=subst(res,t,x,false,contextptr); 1635 return ratnormal(res,contextptr); 1636 } 1637 desolve(const gen & f_orig,const gen & x_orig,const gen & y_orig,int & ordre,vecteur & parameters,GIAC_CONTEXT)1638 gen desolve(const gen & f_orig,const gen & x_orig,const gen & y_orig,int & ordre,vecteur & parameters,GIAC_CONTEXT){ 1639 gen f; 1640 gen x(x_orig),y(y_orig); 1641 if (x.is_symb_of_sommet(at_unquote)) 1642 x=eval(x,1,contextptr); 1643 if (y.is_symb_of_sommet(at_unquote)) 1644 y=eval(y,1,contextptr); 1645 int st=step_infolevel(contextptr); 1646 step_infolevel(0,contextptr); 1647 bool num=false; 1648 gen res=desolve_f(f_orig,x,y,ordre,parameters,f,st,num,contextptr); 1649 if (num) 1650 res=evalf(res,1,contextptr); 1651 step_infolevel(st,contextptr); 1652 return res; 1653 } 1654 1655 // "unary" version _ztrans(const gen & args,GIAC_CONTEXT)1656 gen _ztrans(const gen & args,GIAC_CONTEXT){ 1657 if ( args.type==_STRNG && args.subtype==-1) return args; 1658 if (args.type!=_VECT) 1659 return ztrans(args,vx_var,vx_var,contextptr); 1660 vecteur & v=*args._VECTptr; 1661 int s=int(v.size()); 1662 if (s==2) 1663 return ztrans( v[0],v[1],v[1],contextptr); 1664 if (s!=3) 1665 return gensizeerr(contextptr); 1666 return ztrans( v[0],v[1],v[2],contextptr); 1667 } 1668 static const char _ztrans_s []="ztrans"; 1669 static define_unary_function_eval (__ztrans,&_ztrans,_ztrans_s); 1670 define_unary_function_ptr5( at_ztrans ,alias_at_ztrans,&__ztrans,0,true); 1671 invztranserr(GIAC_CONTEXT)1672 static gen invztranserr(GIAC_CONTEXT){ 1673 return gensizeerr(gettext("Inverse z-transform of non rational functions not implemented or unable to fully factor rational function")); 1674 } 1675 1676 // limited to rational fractions invztrans(const gen & f,const gen & x,const gen & s,GIAC_CONTEXT)1677 gen invztrans(const gen & f,const gen & x,const gen & s,GIAC_CONTEXT){ 1678 if (x.type!=_IDNT) 1679 return gensizeerr(contextptr); 1680 gen t(s); 1681 if (s==x){ 1682 #ifdef GIAC_HAS_STO_38 1683 t=identificateur("s38_"); 1684 #else 1685 t=identificateur(" tinvztrans"); 1686 #endif 1687 } 1688 vecteur varx(lvarx(f,x)); 1689 int varxs=int(varx.size()); 1690 gen res; 1691 if (varxs==0) 1692 res=f*_Kronecker(t,contextptr); 1693 else { 1694 if (varxs>1) 1695 return invztranserr(contextptr); 1696 res=f/x; 1697 vecteur l; 1698 l.push_back(x); // insure x is the main var 1699 l.push_back(t); // s var as second var 1700 l=vecteur(1,l); 1701 alg_lvar(res,l); 1702 vecteur lprime(l); 1703 if (lprime.front().type!=_VECT) return gensizeerr(gettext("desolve.cc/invztrans")); 1704 lprime.front()=cdr_VECT(*(lprime.front()._VECTptr)); 1705 gen glap=e2r(s,l,contextptr); 1706 if (glap.type!=_POLY) return gensizeerr(gettext("desolve.cc/invztrans")); 1707 int dim=int(l.front()._VECTptr->size()); 1708 if (!dim){ 1709 l.erase(l.begin()); 1710 dim=int(l.front()._VECTptr->size()); 1711 } 1712 gen r=e2r(res,l,contextptr); 1713 res=0; 1714 gen r_num,r_den; 1715 fxnd(r,r_num,r_den); 1716 if (r_num.type==_EXT) 1717 return invztranserr(contextptr); 1718 if (r_den.type!=_POLY) 1719 return invztranserr(contextptr); 1720 polynome den(*r_den._POLYptr),num(dim); 1721 if (r_num.type==_POLY) 1722 num=*r_num._POLYptr; 1723 else 1724 num=polynome(r_num,dim); 1725 polynome p_content(lgcd(den)); 1726 den=den/p_content; 1727 factorization vden; gen an; 1728 gen extra_div; 1729 if (!cfactor(den,an,vden,true,extra_div)) 1730 return invztranserr(contextptr); 1731 vector< pf<gen> > pfde_VECT; 1732 polynome ipnum(dim),ipden(dim); 1733 partfrac(num,den,vden,pfde_VECT,ipnum,ipden); 1734 if (!is_zero(ipnum)) 1735 *logptr(contextptr) << gettext("Warning, z*argument has a non-zero integral part") << '\n'; 1736 vector< pf<gen> >::iterator it=pfde_VECT.begin(); 1737 vector< pf<gen> >::const_iterator itend=pfde_VECT.end(); 1738 gen a,A,B; 1739 polynome b,c; 1740 for (;it!=itend;++it){ 1741 if (it->fact.lexsorted_degree()>1) 1742 return invztranserr(contextptr); 1743 findde(it->fact,b,c); 1744 a=-gen(c)/gen(b); // pole 1745 B=r2e(Tfirstcoeff(it->den),l,contextptr); 1746 if (is_zero(a)){ 1747 int mult=it->mult; 1748 gen res0; 1749 vecteur vnum; 1750 polynome2poly1(it->num,1,vnum); 1751 for (int i=0;i<mult;++i){ 1752 res0 += r2e(vnum[i],lprime,contextptr)*symbolic(at_Kronecker,s-i); // symb_when(symb_equal(s,i),1,0) will not be handled correctly by ztrans 1753 } 1754 res += res0/B; 1755 } 1756 else { 1757 // it->num/it->den in terms of 1/(z-a), a/(z-a)^2, a^2/(z-a)^3, etc. 1758 gen cur=r2e(it->num,l,contextptr); 1759 A=r2e(a,lprime,contextptr); 1760 gen z_minus_a=x-A,res0; 1761 for (int i=it->mult-1;i>=0;--i){ 1762 gen tmp=_quorem(makesequence(cur,z_minus_a,x),contextptr); 1763 if (is_undef(tmp)) return tmp; 1764 gen rem=tmp[1]; 1765 cur=tmp[0]; 1766 rem=rem/pow(A,i,contextptr)/factorial(i); 1767 for (int j=0;j<i;++j) 1768 rem = rem * (s-j); 1769 res0 += rem; 1770 } 1771 res0 = res0 * pow(A,s,contextptr); 1772 res += res0/B; 1773 } 1774 } 1775 res=res/r2e(p_content,l,contextptr); 1776 } 1777 if (s==x) 1778 res=subst(res,t,x,false,contextptr); 1779 res=ratnormal(res,contextptr); 1780 // replace discrete Kronecker by Heaviside in some very simple situations 1781 vecteur vD=lop(res,at_Kronecker); 1782 gen A,B,a,b; 1783 if (vD.size()==1 && is_linear_wrt(res,vD.front(),A,B,contextptr) && is_linear_wrt(vD.front()._SYMBptr->feuille,s,a,b,contextptr)){ 1784 // res==A*Kronecker(a*x+b)+B 1785 if (is_one(a) && is_zero(b)){ 1786 gen B0=subst(B,s,0,false,contextptr); 1787 if (is_zero(ratnormal(B0+A,contextptr))) 1788 res=B*symbolic(at_Heaviside,s-1); 1789 } 1790 } 1791 return res; 1792 } 1793 _invztrans(const gen & args,GIAC_CONTEXT)1794 gen _invztrans(const gen & args,GIAC_CONTEXT){ 1795 if ( args.type==_STRNG && args.subtype==-1) return args; 1796 if (args.type!=_VECT) 1797 return invztrans(args,vx_var,vx_var,contextptr); 1798 vecteur & v=*args._VECTptr; 1799 int s=int(v.size()); 1800 if (s==2) 1801 return invztrans( v[0],v[1],v[1],contextptr); 1802 if (s!=3) 1803 return gensizeerr(contextptr); 1804 return invztrans( v[0],v[1],v[2],contextptr); 1805 } 1806 static const char _invztrans_s []="invztrans"; 1807 static define_unary_function_eval (__invztrans,&_invztrans,_invztrans_s); 1808 define_unary_function_ptr5( at_invztrans ,alias_at_invztrans,&__invztrans,0,true); 1809 _Kronecker(const gen & args,GIAC_CONTEXT)1810 gen _Kronecker(const gen & args,GIAC_CONTEXT){ 1811 if ( args.type==_STRNG && args.subtype==-1) return args; 1812 if (args.type==_VECT) 1813 return apply(args,_Kronecker,contextptr); 1814 if (!is_integer(args)) 1815 return symbolic(at_Kronecker,args); 1816 if (is_zero(args)) 1817 return 1; 1818 else 1819 return 0; 1820 } 1821 static const char _Kronecker_s []="Kronecker"; 1822 static define_unary_function_eval (__Kronecker,&_Kronecker,_Kronecker_s); 1823 define_unary_function_ptr5( at_Kronecker ,alias_at_Kronecker,&__Kronecker,0,true); 1824 1825 #ifndef USE_GMP_REPLACEMENTS 1826 // this code slice (c) 2018 Luka Marohnić 1827 /* returns the coefficient of t in (fraction, Taylor or Laurent) expansion e */ expansion_coeff(const gen & e,const gen & t,GIAC_CONTEXT)1828 gen expansion_coeff(const gen &e,const gen &t,GIAC_CONTEXT) { 1829 gen ret(0),g; 1830 if (e.is_symb_of_sommet(at_plus) && e._SYMBptr->feuille.type==_VECT) { 1831 const vecteur &feu=*e._SYMBptr->feuille._VECTptr; 1832 for (const_iterateur it=feu.begin();it!=feu.end();++it) { 1833 g=_ratnormal(*it/t,contextptr); 1834 if (_evalf(g,contextptr).type==_DOUBLE_) { 1835 ret=g; 1836 break; 1837 } 1838 } 1839 } else { 1840 g=_ratnormal(e/t,contextptr); 1841 if (_evalf(g,contextptr).type==_DOUBLE_) 1842 ret=g; 1843 } 1844 return ret; 1845 } 1846 kovacic_iscase1(const vecteur & poles,int dinf)1847 bool kovacic_iscase1(const vecteur &poles,int dinf) { 1848 if (dinf%2!=0 && dinf<=2) 1849 return false; 1850 for (const_iterateur it=poles.begin();it!=poles.end();++it) { 1851 int order=it->_VECTptr->back().val; 1852 if (order%2!=0 && order!=1) 1853 return false; 1854 } 1855 return true; 1856 } 1857 kovacic_iscase2(const vecteur & poles)1858 bool kovacic_iscase2(const vecteur &poles) { 1859 for (const_iterateur it=poles.begin();it!=poles.end();++it) { 1860 int order=it->_VECTptr->back().val; 1861 if (order==2 || (order>2 && order%2!=0)) 1862 return true; 1863 } 1864 return false; 1865 } 1866 kovacic_iscase3(const gen & cpfr,const gen & x,const vecteur & poles,int dinf,GIAC_CONTEXT)1867 bool kovacic_iscase3(const gen &cpfr,const gen &x,const vecteur &poles,int dinf,GIAC_CONTEXT) { 1868 if (dinf<2) 1869 return false; 1870 vecteur alpha,beta; 1871 gen a,b,g(0),p; 1872 int d; 1873 for (const_iterateur it=poles.begin();it!=poles.end();++it) { 1874 p=it->_VECTptr->front(); 1875 d=it->_VECTptr->back().val; 1876 if (d>2) 1877 return false; 1878 if (d>1) { 1879 alpha.push_back(a=expansion_coeff(cpfr,pow(x-p,-2),contextptr)); 1880 g+=a; 1881 a=_eval(sqrt(1+4*a,contextptr),contextptr); 1882 if (!_numer(a,contextptr).is_integer() || !_denom(a,contextptr).is_integer()) 1883 return false; 1884 } 1885 beta.push_back(b=expansion_coeff(cpfr,_inv(x-p,contextptr),contextptr)); 1886 g+=b*p; 1887 } 1888 g=_eval(sqrt(1+4*g,contextptr),contextptr); 1889 return is_zero(_ratnormal(_sum(beta,contextptr),contextptr)) && 1890 _numer(g,contextptr).is_integer() && _denom(g,contextptr).is_integer(); 1891 } 1892 build_E_families(const gen_map & E,const vecteur & cv,vecteur & family,matrice & families)1893 void build_E_families(const gen_map &E,const vecteur &cv,vecteur &family,matrice &families) { 1894 int i=family.size(); 1895 if (i>=int(cv.size())) 1896 return; 1897 const vecteur ev=*E.find(cv[i])->second._VECTptr; 1898 for (const_iterateur it=ev.begin();it!=ev.end();++it) { 1899 family.push_back(*it); 1900 if (family.size()==cv.size()) 1901 families.push_back(family); 1902 else build_E_families(E,cv,family,families); 1903 family.pop_back(); 1904 } 1905 } 1906 create_identifiers(vecteur & vars,int n)1907 void create_identifiers(vecteur &vars,int n) { 1908 vars.reserve(n); 1909 stringstream ss; 1910 for (int i=0;i<n;++i) { 1911 ss.str(""); 1912 ss << " cf" << i; 1913 vars.push_back(identificateur(ss.str().c_str())); 1914 } 1915 } 1916 strip_abs(const gen & g)1917 gen strip_abs(const gen &g) { 1918 if (g.is_symb_of_sommet(at_abs)) 1919 return g._SYMBptr->feuille; 1920 if (g.type==_SYMB) { 1921 gen args; 1922 if (g._SYMBptr->feuille.type==_VECT) { 1923 args=vecteur(0); 1924 const vecteur &feu=*g._SYMBptr->feuille._VECTptr; 1925 for (const_iterateur it=feu.begin();it!=feu.end();++it) { 1926 args._VECTptr->push_back(strip_abs(*it)); 1927 } 1928 } else args=strip_abs(g._SYMBptr->feuille); 1929 return symbolic(g._SYMBptr->sommet,args); 1930 } 1931 return g; 1932 } 1933 explnsimp(const gen & g,GIAC_CONTEXT)1934 gen explnsimp(const gen &g,GIAC_CONTEXT) { 1935 gen e=expand(strip_abs(g),contextptr); 1936 e=symbolic(at_exp2pow,symbolic(at_expexpand,symbolic(at_lncollect,e))); 1937 return ratnormal(_lin(_eval(e,contextptr),contextptr),contextptr); 1938 } 1939 isroot(const gen & g,gen & deg,GIAC_CONTEXT)1940 bool isroot(const gen &g,gen °,GIAC_CONTEXT) { 1941 if (!g.is_symb_of_sommet(at_pow)) 1942 return false; 1943 const gen &pw=g._SYMBptr->feuille._VECTptr->at(1); 1944 return (deg=_inv(pw,contextptr)).is_integer() && !is_minus_one(deg); 1945 } 1946 partialrad(const gen & g,const gen & deg,gen & outside,gen & inside,bool isdenom,GIAC_CONTEXT)1947 void partialrad(const gen &g,const gen °,gen &outside,gen &inside,bool isdenom,GIAC_CONTEXT) { 1948 gen s; 1949 if (g.is_integer() && (s=_pow(makesequence(g,_inv(deg,contextptr)),contextptr)).is_integer()) { 1950 outside=outside*(isdenom?_inv(s,contextptr):s); 1951 } else if (g.is_symb_of_sommet(at_prod) && g._SYMBptr->feuille.type==_VECT) { 1952 vecteur &fv=*g._SYMBptr->feuille._VECTptr,qr; 1953 for (const_iterateur it=fv.begin();it!=fv.end();++it) { 1954 if (it->is_integer()) { 1955 s=_pow(makesequence(*it,_inv(deg,contextptr)),contextptr); 1956 outside=outside*(isdenom?_inv(s,contextptr):s); 1957 continue; 1958 } else if (it->is_symb_of_sommet(at_pow)) { 1959 const gen &pw=it->_SYMBptr->feuille._VECTptr->at(1); 1960 const gen &b=it->_SYMBptr->feuille._VECTptr->front(); 1961 if (pw.is_integer()) { 1962 qr=*_iquorem(makesequence(pw,deg),contextptr)._VECTptr; 1963 if (isdenom) { 1964 outside=outside/_pow(makesequence(b,qr.front()),contextptr); 1965 inside=inside/_pow(makesequence(b,qr.back()),contextptr); 1966 } else { 1967 outside=outside*_pow(makesequence(b,qr.front()),contextptr); 1968 inside=inside*_pow(makesequence(b,qr.back()),contextptr); 1969 } 1970 continue; 1971 } 1972 } else if (it->is_symb_of_sommet(at_inv)) { 1973 partialrad(it->_SYMBptr->feuille,deg,outside,inside,!isdenom,contextptr); 1974 continue; 1975 } 1976 inside=inside*(isdenom?_inv(*it,contextptr):*it); 1977 } 1978 } else if (g.is_symb_of_sommet(at_inv)) 1979 partialrad(g._SYMBptr->feuille,deg,outside,inside,!isdenom,contextptr); 1980 else inside=inside*(isdenom?_inv(g,contextptr):g); 1981 } 1982 radsimp(const gen & g,GIAC_CONTEXT)1983 gen radsimp(const gen &g,GIAC_CONTEXT) { 1984 gen deg; 1985 if (g.type==_VECT) { 1986 vecteur ret; 1987 for (const_iterateur it=g._VECTptr->begin();it!=g._VECTptr->end();++it) { 1988 ret.push_back(radsimp(*it,contextptr)); 1989 } 1990 return change_subtype(ret,g.subtype); 1991 } 1992 if (isroot(g,deg,contextptr)) { 1993 gen radic=_collect(radsimp(g._SYMBptr->feuille._VECTptr->front(),contextptr),contextptr); 1994 gen inside(1),outside(1),inum; 1995 partialrad(radic,deg,outside,inside,false,contextptr); 1996 gen ideg=_inv(deg,contextptr); 1997 inside=_eval(symb_normal(inside),contextptr); 1998 if (!(inum=_eval(_pow(makesequence(_numer(inside,contextptr),ideg),contextptr),contextptr)).is_symb_of_sommet(at_pow) || 1999 inum._SYMBptr->feuille._VECTptr->at(1).is_integer()) 2000 return _collect(outside,contextptr)*inum/_pow(makesequence(_collect(_denom(inside,contextptr),contextptr),ideg),contextptr); 2001 return _collect(outside,contextptr)*_pow(makesequence(_collect(inside,contextptr),ideg),contextptr); 2002 } 2003 if (g.is_symb_of_sommet(at_prod) && g._SYMBptr->feuille.type==_VECT) { 2004 vecteur &feu=*g._SYMBptr->feuille._VECTptr,degv; 2005 gen den(1); 2006 for (const_iterateur it=feu.begin();it!=feu.end();++it) { 2007 if (it->is_symb_of_sommet(at_pow)) { 2008 gen dg=_inv(it->_SYMBptr->feuille._VECTptr->at(1),contextptr); 2009 if (is_greater(dg,2,contextptr)) 2010 degv.push_back(dg); 2011 } else if (it->is_symb_of_sommet(at_inv)) 2012 den=den*it->_SYMBptr->feuille; 2013 } 2014 if (den.is_symb_of_sommet(at_prod) && den._SYMBptr->feuille.type==_VECT) { 2015 vecteur &dfeu=*den._SYMBptr->feuille._VECTptr; 2016 for (const_iterateur it=dfeu.begin();it!=dfeu.end();++it) { 2017 if (it->is_symb_of_sommet(at_pow)) { 2018 gen dg=_inv(it->_SYMBptr->feuille._VECTptr->at(1),contextptr); 2019 if (is_greater(dg,2,contextptr)) 2020 degv.push_back(dg); 2021 } 2022 } 2023 } 2024 gen gd=_lcm(degv,contextptr); 2025 gen p=_collect(ratnormal(pow(g,gd.val),contextptr),contextptr); 2026 if (p.is_symb_of_sommet(at_pow)) { 2027 gen ppw=p._SYMBptr->feuille._VECTptr->at(1); 2028 if (ppw.is_integer()) 2029 gd=_eval(gd/ppw,contextptr); 2030 } 2031 gen ret=_pow(makesequence(p,_inv(gd,contextptr)),contextptr); 2032 return isroot(ret,deg,contextptr)?radsimp(ret,contextptr):ratnormal(ret,contextptr); 2033 } 2034 if (g.type==_SYMB) { 2035 gen &feu=g._SYMBptr->feuille; 2036 if (feu.type==_VECT) { 2037 vecteur res; 2038 for (const_iterateur it=feu._VECTptr->begin();it!=feu._VECTptr->end();++it) { 2039 res.push_back(radsimp(*it,contextptr)); 2040 } 2041 return symbolic(g._SYMBptr->sommet,change_subtype(res,feu.subtype)); 2042 } 2043 return symbolic(g._SYMBptr->sommet,radsimp(feu,contextptr)); 2044 } 2045 return g; 2046 } 2047 strip_gcd(const vecteur & v,GIAC_CONTEXT)2048 vecteur strip_gcd(const vecteur &v,GIAC_CONTEXT) { 2049 gen g1=_gcd(_apply(makesequence(at_numer,v),contextptr),contextptr); 2050 gen g2=_gcd(_apply(makesequence(at_denom,v),contextptr),contextptr); 2051 return *_collect(_ratnormal(multvecteur(g2/g1,v),contextptr),contextptr)._VECTptr; 2052 } 2053 ratsimp_nonexp(const gen & g,GIAC_CONTEXT)2054 gen ratsimp_nonexp(const gen &g,GIAC_CONTEXT) { 2055 if (g.type==_VECT) { 2056 vecteur res; 2057 for (const_iterateur it=g._VECTptr->begin();it!=g._VECTptr->end();++it) { 2058 res.push_back(ratsimp_nonexp(*it,contextptr)); 2059 } 2060 return change_subtype(res,g.subtype); 2061 } 2062 if (g.is_symb_of_sommet(at_plus) && g._SYMBptr->feuille==_VECT) { 2063 vecteur &terms=*g._SYMBptr->feuille._VECTptr; 2064 gen res(0); 2065 for (const_iterateur it=terms.begin();it!=terms.end();++it) { 2066 res+=ratsimp_nonexp(*it,contextptr); 2067 } 2068 return res; 2069 } 2070 if (g.is_symb_of_sommet(at_prod) && g._SYMBptr->feuille.type==_VECT) { 2071 vecteur &facs=*g._SYMBptr->feuille._VECTptr; 2072 gen e(1),ne(1); 2073 for (const_iterateur it=facs.begin();it!=facs.end();++it) { 2074 if (it->is_symb_of_sommet(at_exp)) 2075 e=*it*e; 2076 else ne=*it*ne; 2077 } 2078 return _ratnormal(ne,contextptr)*e; 2079 } 2080 return g; 2081 } 2082 fullsimp(const gen & g,GIAC_CONTEXT)2083 gen fullsimp(const gen &g,GIAC_CONTEXT) { 2084 return ratsimp_nonexp(_collect(radsimp(explnsimp(exp(_ratnormal(g,contextptr),contextptr), 2085 contextptr),contextptr),contextptr),contextptr); 2086 } 2087 2088 /* 2089 * This routine solves the general homogeneous linear second-order ODE 2090 * y''=r(t)*y, where r is a non-constant rational function, using Kovacic's 2091 * algorithm (https://core.ac.uk/download/pdf/82509765.pdf). A list of 2092 * solutions is returned (possibly empty). 2093 */ kovacicsols(const gen & r_orig,const gen & x,const gen & dy_coeff,GIAC_CONTEXT)2094 gen kovacicsols(const gen &r_orig,const gen &x,const gen &dy_coeff,GIAC_CONTEXT) { 2095 gen r=_ratnormal(r_orig,contextptr),inf=symbolic(at_plus,_IDNT_infinity()); 2096 gen s=_numer(r,contextptr),t=_denom(r,contextptr),a,b,c,e,w=identificateur("omega_"); 2097 int ds=_degree(makesequence(s,x),contextptr).val; 2098 int dt=_degree(makesequence(t,x),contextptr).val,dinf=dt-ds,order,nu; 2099 vecteur poles=*_roots(makesequence(t,x),contextptr)._VECTptr,solutions(0); 2100 gen cpfr=_cpartfrac(makesequence(r,x),contextptr); 2101 gen laur=_series(makesequence(r,x,inf,1,_POLY1__VECT),contextptr); 2102 bool success=false; 2103 if (kovacic_iscase1(poles,dinf)) { 2104 cerr << "Case 1 of Kovacic algorithm" << '\n'; 2105 /* step 1 */ 2106 gen_map alpha_plus,alpha_minus,sqrt_r; 2107 gen alpha_inf_plus,alpha_inf_minus; 2108 for (const_iterateur it=poles.begin();it!=poles.end();++it) { 2109 c=it->_VECTptr->front(); 2110 order=it->_VECTptr->back().val; 2111 if (order==1) { 2112 sqrt_r[c]=0; 2113 alpha_plus[c]=alpha_minus[c]=1; 2114 } else if (order==2) { 2115 sqrt_r[c]=0; 2116 b=expansion_coeff(cpfr,pow(x-c,-2),contextptr); 2117 alpha_plus[c]=(1+sqrt(1+4*b,contextptr))/2; 2118 alpha_minus[c]=(1-sqrt(1+4*b,contextptr))/2; 2119 } else if (order%2==0 && order>=4) { 2120 nu=order/2; 2121 e=_series(makesequence(sqrt(r,contextptr),x,c,1,_POLY1__VECT),contextptr); 2122 for (int i=2;i<=nu;++i) { 2123 gen cf=expansion_coeff(e,pow(x-c,-i),contextptr); 2124 if (i==nu) 2125 a=cf; 2126 sqrt_r[c]+=cf/pow(x-c,i); 2127 } 2128 if (is_zero(a)) 2129 return false; 2130 b=expansion_coeff(cpfr,pow(x-c,-nu-1),contextptr); 2131 b-=expansion_coeff(_cpartfrac(makesequence(sq(sqrt_r[c]),x),contextptr),pow(x-c,-nu-1),contextptr); 2132 alpha_plus[c]=(nu+b/a)/2; 2133 alpha_minus[c]=(nu-b/a)/2; 2134 } else assert(false); 2135 } 2136 if (dinf>2) { 2137 sqrt_r[inf]=0; 2138 alpha_inf_plus=0; 2139 alpha_inf_minus=1; 2140 } else if (dinf==2) { 2141 sqrt_r[inf]=0; 2142 b=_lcoeff(makesequence(s,x),contextptr)/_lcoeff(makesequence(t,x),contextptr); 2143 alpha_inf_plus=(1+sqrt(1+4*b,contextptr))/2; 2144 alpha_inf_minus=(1-sqrt(1+4*b,contextptr))/2; 2145 } else if (dinf%2==0 && dinf<=0) { 2146 nu=-dinf/2; 2147 e=_series(makesequence(sqrt(r,contextptr),x,inf,1,_POLY1__VECT),contextptr); 2148 for (int i=0;i<=nu;++i) { 2149 gen cf=expansion_coeff(e,pow(x,i),contextptr); 2150 if (i==nu) 2151 a=cf; 2152 sqrt_r[inf]+=cf*pow(x,i); 2153 } 2154 if (is_zero(a)) 2155 return false; 2156 b=expansion_coeff(_propfrac(makesequence(r,x),contextptr),pow(x,nu-1),contextptr); 2157 b-=expansion_coeff(expand(sq(sqrt_r[inf]),contextptr),pow(x,nu-1),contextptr); 2158 alpha_inf_plus=(b/a-nu)/2; 2159 alpha_inf_minus=(-b/a-nu)/2; 2160 } else assert(false); 2161 /* step 2 */ 2162 int np=poles.size()+1,N=(1<<np),sc,d,maxd=0; // N=std::pow(2.0,np) 2163 vecteur fam,cv,vars,v; 2164 gen fd,fw,P; 2165 for (gen_map::const_iterator it=alpha_plus.begin();it!=alpha_plus.end();++it) { 2166 cv.push_back(it->first); 2167 } 2168 for (int i=0;i<N;++i) { 2169 fd=(i & 1)!=0?alpha_inf_plus:alpha_inf_minus; 2170 for (int j=1;j<np;++j) { 2171 fd-=(i & (1<<j))!=0?alpha_plus[cv[j-1]]:alpha_minus[cv[j-1]]; // i & (int)std::pow(2.0,j) 2172 } 2173 fd=_ratnormal(fd,contextptr); 2174 if (fd.is_integer() && is_positive(fd,contextptr)) { 2175 fw=(i & 1)!=0?sqrt_r[inf]:-sqrt_r[inf]; 2176 for (int j=1;j<np;++j) { 2177 sc=(i & (1<<j)); // i & (int)std::pow(2.0,j) 2178 fw+=sc!=0?sqrt_r[cv[j-1]]:-sqrt_r[cv[j-1]]; 2179 fw+=(sc!=0?alpha_plus[cv[j-1]]:alpha_minus[cv[j-1]])/(x-cv[j-1]); 2180 } 2181 fam.push_back(makevecteur(fd,fw)); 2182 maxd=std::max(maxd,fd.val); 2183 } 2184 } 2185 /* step 3 */ 2186 create_identifiers(vars,maxd); 2187 for (const_iterateur it=fam.begin();it!=fam.end() && !success;++it) { 2188 d=it->_VECTptr->front().val; fw=it->_VECTptr->back(); 2189 v=vecteur(vars.begin(),vars.begin()+d); 2190 P=0; 2191 for (int i=0;i<d;++i) P+=v[i]*pow(x,i); 2192 P+=pow(x,d); 2193 e=_numer(_derive(makesequence(P,x,2),contextptr)+2*fw*_derive(makesequence(P,x),contextptr)+ 2194 (_derive(makesequence(fw,x),contextptr)+sq(fw)-r)*P,contextptr); 2195 gen lsol=_solve(makesequence(_coeff(makesequence(e,x),contextptr),v),contextptr); 2196 if (lsol.type==_VECT && !lsol._VECTptr->empty()) { 2197 lsol=_subst(makesequence(lsol._VECTptr->front(),v,vecteur(d,0)),contextptr); 2198 solutions.push_back(_subst(makesequence(P,v,lsol),contextptr)* 2199 fullsimp(_int(makesequence(fw-dy_coeff/2,x),contextptr),contextptr)); 2200 success=true; 2201 } 2202 } 2203 } 2204 if (!success && kovacic_iscase2(poles)) { 2205 cerr << "Case 2 of Kovacic algorithm" << '\n'; 2206 /* step 1 */ 2207 gen_map E; 2208 for (const_iterateur it=poles.begin();it!=poles.end();++it) { 2209 c=it->_VECTptr->front(); 2210 order=it->_VECTptr->back().val; 2211 if (order==1) 2212 E[c]=vecteur(1,4); 2213 else if (order==2) { 2214 b=expansion_coeff(cpfr,pow(x-c,-2),contextptr); 2215 E[c]=vecteur(0); 2216 for (int k=-1;k<=1;++k) { 2217 gen tmp=_eval(2+2*k*sqrt(1+4*b,contextptr),contextptr); 2218 if (tmp.is_integer()) 2219 E[c]._VECTptr->push_back(tmp); 2220 } 2221 } else if (order>2) 2222 E[c]=vecteur(1,order); 2223 } 2224 vecteur Einf(0); 2225 if (dinf>2) 2226 Einf=makevecteur(0,2,4); 2227 else if (dinf==2) { 2228 b=expansion_coeff(laur,pow(x,-2),contextptr); 2229 for (int k=-1;k<=1;++k) { 2230 gen tmp=_eval(2+2*k*sqrt(1+4*b,contextptr),contextptr); 2231 if (tmp.is_integer()) 2232 Einf.push_back(tmp); 2233 } 2234 } else if (dinf<2) 2235 Einf.push_back(dinf); 2236 /* step 2 */ 2237 vecteur family,families,fam,cv,vars,v; 2238 for (gen_map::const_iterator it=E.begin();it!=E.end();++it) { 2239 cv.push_back(it->first); 2240 } 2241 build_E_families(E,cv,family,families); 2242 int maxdeg=0,deg; 2243 for (const_iterateur it=Einf.begin();it!=Einf.end();++it) { 2244 if (families.empty()) { 2245 e=*it/2; 2246 if (e.is_integer() && is_positive(e,contextptr)) { 2247 fam.push_back(makevecteur(e,0)); 2248 maxdeg=std::max(maxdeg,e.val); 2249 } 2250 } 2251 for (const_iterateur jt=families.begin();jt!=families.end();++jt) { 2252 gen th(0); 2253 bool discard=is_one(_even(*it,contextptr)); 2254 const vecteur &fm=*(jt->_VECTptr); 2255 for (const_iterateur kt=fm.begin();kt!=fm.end();++kt) { 2256 th+=(*kt)/(x-cv[kt-fm.begin()]); 2257 if (is_zero(_even(*kt,contextptr))) 2258 discard=false; 2259 } 2260 //if (discard) continue; 2261 e=_eval(*it-_sum(fm,contextptr),contextptr)/2; 2262 if (e.is_integer() && is_positive(e,contextptr)) { 2263 fam.push_back(makevecteur(e,th/2)); 2264 maxdeg=std::max(maxdeg,e.val); 2265 } 2266 } 2267 } 2268 /* step 3 */ 2269 create_identifiers(vars,maxdeg); 2270 gen P,th,dth; 2271 for (const_iterateur it=fam.begin();it!=fam.end() && !success;++it) { 2272 deg=it->_VECTptr->front().val; th=it->_VECTptr->back(); 2273 v=vecteur(vars.begin(),vars.begin()+deg); 2274 P=0; 2275 for (int i=0;i<deg;++i) P+=v[i]*pow(x,i); 2276 P+=pow(x,deg); 2277 dth=_derive(makesequence(th,x),contextptr); 2278 e=_derive(makesequence(P,x,3),contextptr)+3*th*_derive(makesequence(P,x,2),contextptr)+ 2279 (3*sq(th)+3*dth-4*r)*_derive(makesequence(P,x),contextptr)+ 2280 (_derive(makesequence(th,x,2),contextptr)+3*th*dth+pow(th,3)-4*r*th-2*_derive(makesequence(r,x),contextptr))*P; 2281 e=_numer(e,contextptr); 2282 gen cfs=deg==0?undef:_solve(makesequence(_coeff(makesequence(e,x),contextptr),v),contextptr); 2283 if (deg==0?is_zero(e):(cfs.type==_VECT && !cfs._VECTptr->empty())) { 2284 if (deg>0) { 2285 cfs=_subst(makesequence(cfs._VECTptr->front(),v,vecteur(deg,0)),contextptr); 2286 P=_subst(makesequence(P,v,cfs),contextptr); 2287 } 2288 gen ph=th+_derive(makesequence(P,x),contextptr)/P; 2289 gen qsol=_solve(makesequence(symb_equal(sq(w)-w*ph+(_derive(makesequence(ph,x),contextptr)/2+sq(ph)/2-r),0),w),contextptr); 2290 if (qsol.type==_VECT) for (const_iterateur jt=qsol._VECTptr->begin();jt!=qsol._VECTptr->end();++jt) { 2291 solutions.push_back(fullsimp(_int(makesequence(*jt-dy_coeff/2,x),contextptr),contextptr)); 2292 success=true; 2293 } 2294 } 2295 } 2296 } 2297 if (!success && kovacic_iscase3(cpfr,x,poles,dinf,contextptr)) { 2298 cerr << "Case 3 of Kovacic algorithm" << '\n'; 2299 vector<int> nv=vecteur_2_vector_int(makevecteur(4,6,12)); 2300 for (vector<int>::const_iterator nt=nv.begin();nt!=nv.end();++nt) { 2301 int n=*nt; 2302 /* step 1 */ 2303 gen_map E; 2304 for (const_iterateur it=poles.begin();it!=poles.end();++it) { 2305 c=it->_VECTptr->front(); 2306 order=it->_VECTptr->back().val; 2307 if (order==1) 2308 E[c]=vecteur(1,12); 2309 else if (order==2) { 2310 a=expansion_coeff(cpfr,pow(x-c,-2),contextptr); 2311 E[c]=vecteur(0); 2312 for (int k=-n/2;k<=n/2;++k) { 2313 gen tmp=_eval(6+k*(12/n)*sqrt(1+4*a,contextptr),contextptr); 2314 if (tmp.is_integer()) 2315 E[c]._VECTptr->push_back(tmp); 2316 } 2317 } 2318 } 2319 vecteur Einf(0); 2320 b=dinf>2?gen(0):expansion_coeff(laur,pow(x,-2),contextptr); 2321 for (int k=-n/2;k<=n/2;++k) { 2322 gen tmp=_eval(6+k*(12/n)*sqrt(1+4*b,contextptr),contextptr); 2323 if (tmp.is_integer()) 2324 Einf.push_back(tmp); 2325 } 2326 /* step 2 */ 2327 vecteur family,families,fam,cv,vars,v; 2328 for (gen_map::const_iterator it=E.begin();it!=E.end();++it) { 2329 cv.push_back(it->first); 2330 } 2331 build_E_families(E,cv,family,families); 2332 int maxdeg=0,deg; 2333 for (const_iterateur it=Einf.begin();it!=Einf.end();++it) { 2334 if (families.empty()) { 2335 e=*it*gen(n)/12; 2336 if (e.is_integer() && is_positive(e,contextptr)) { 2337 fam.push_back(makevecteur(e,0)); 2338 maxdeg=std::max(maxdeg,e.val); 2339 } 2340 } 2341 for (const_iterateur jt=families.begin();jt!=families.end();++jt) { 2342 gen th(0); 2343 const vecteur &fm=*(jt->_VECTptr); 2344 for (const_iterateur kt=fm.begin();kt!=fm.end();++kt) { 2345 th+=(*kt)/(x-cv[kt-fm.begin()]); 2346 } 2347 e=gen(n)*_eval(*it-_sum(*jt,contextptr),contextptr)/12; 2348 if (e.is_integer() && is_positive(e,contextptr)) { 2349 fam.push_back(makevecteur(e,gen(n)*th/12)); 2350 maxdeg=std::max(maxdeg,e.val); 2351 } 2352 } 2353 } 2354 /* step 3 */ 2355 gen S(1); 2356 for (const_iterateur it=cv.begin();it!=cv.end();++it) { 2357 S=S*(x-*it); 2358 } 2359 gen dS=_derive(makesequence(S,x),contextptr),th; 2360 vecteur P(n+2,0); 2361 create_identifiers(vars,maxdeg); 2362 for (const_iterateur it=fam.begin();it!=fam.end();++it) { 2363 deg=it->_VECTptr->front().val; th=it->_VECTptr->back(); 2364 v=vecteur(vars.begin(),vars.begin()+deg); 2365 for (int i=0;i<deg;++i) P[n+1]-=v[i]*pow(x,i); 2366 P[n+1]-=pow(x,deg); 2367 P[n]=-S*(_derive(makesequence(P[n+1],x),contextptr)+th*P[n+1]),contextptr; 2368 if (P[n].type==_SYMB) P[n]=_collect(P[n],contextptr); 2369 for (int i=n;i-->0;) { 2370 P[i]=-S*_derive(makesequence(P[i+1],x),contextptr)+((n-i)*dS-S*th)*P[i+1]-(n-i)*(i+1)*sq(S)*r*P[i+2]; 2371 if (P[i].type==_SYMB) P[i]=_collect(P[i],contextptr); 2372 } 2373 gen cfs; 2374 if ((deg==0 && is_zero(P[0])) || 2375 (deg>0 && (cfs=_solve(makesequence(_coeff(makesequence(P[0],x),contextptr),v),contextptr)).type==_VECT && 2376 !cfs._VECTptr->empty())) { 2377 if (deg>0) { 2378 cfs=_subst(makesequence(cfs._VECTptr->front(),v,vecteur(deg,0)),contextptr); 2379 P=*_subst(makesequence(P,v,cfs),contextptr)._VECTptr; 2380 } 2381 vecteur ac(n+1); 2382 for (int i=0;i<=n;++i) { 2383 ac[i]=_collect(pow(S,i)*P[i+1]/_factorial(n-i,contextptr),contextptr); 2384 } 2385 *logptr(contextptr) << "Warning: outputting the algebraic expression for ω" << '\n'; 2386 ac=strip_gcd(ac,contextptr); 2387 gen omg=pow(w,4)*ac[4]+pow(w,3)*ac[3]+pow(w,2)*ac[2]+w*ac[1]+ac[0]; 2388 if (!is_zero(dy_coeff)) { 2389 vecteur C=*_coeff(makesequence(_subst(makesequence(omg,w,w+dy_coeff/2),contextptr),w),contextptr)._VECTptr; 2390 for (int i=0;i<=n;++i) { 2391 ac[i]=_collect(_ratnormal(C[i],contextptr),contextptr); 2392 } 2393 ac=strip_gcd(ac,contextptr); 2394 omg=pow(w,4)*ac[4]+pow(w,3)*ac[3]+pow(w,2)*ac[2]+w*ac[1]+ac[0]; 2395 } 2396 return omg; 2397 } 2398 } 2399 } 2400 } 2401 return solutions; 2402 } 2403 2404 /* 2405 * Return the solution(s) of a second-order linear homogeneous ODE using 2406 * Kovacic's algorithm. The first argument is the ODE a(x)*y''+b(x)*y'+c(x)*y=0 2407 * itself, which may be given as an expression (left-hand side), an equation or 2408 * a list [a,b,c]. The functions a, b and c must be rational in x. The second 2409 * and third (both optional) arguments are the independent variable x and the 2410 * dependent variable y, respectively. By default, the symbols "x" and "y" are 2411 * used. 2412 */ _kovacicsols(const gen & g,GIAC_CONTEXT)2413 gen _kovacicsols(const gen &g,GIAC_CONTEXT) { 2414 if (g.type==_STRNG && g.subtype==-1) return g; 2415 gen x=identificateur("x"),y=identificateur("y"),eq,p(0),q(0),r(0); 2416 if (g.type!=_VECT || g.subtype!=_SEQ__VECT) { 2417 eq=g; 2418 } else if (g.subtype==_SEQ__VECT) { 2419 vecteur &gv=*g._VECTptr; 2420 if (gv.size()<2) return gensizeerr(contextptr); 2421 eq=gv.front(); 2422 if (gv.size()==2) { 2423 if (eq.type==_VECT && gv.back().type==_IDNT) 2424 x=gv.back(); 2425 else if (gv.back().is_symb_of_sommet(at_of)) { 2426 y=gv.back()._SYMBptr->feuille._VECTptr->front(); 2427 x=gv.back()._SYMBptr->feuille._VECTptr->back(); 2428 } else return gensizeerr(contextptr); 2429 } else { 2430 if (eq.type==_VECT) 2431 return gensizeerr(contextptr); 2432 x=gv[1]; 2433 y=gv[2]; 2434 } 2435 if (y.type!=_IDNT || x.type!=_IDNT) 2436 return gensizeerr(contextptr); 2437 } 2438 eq=idnteval(eq,contextptr); 2439 if (eq.type==_VECT) { 2440 vecteur &cfs=*eq._VECTptr; 2441 if (cfs.size()!=3) 2442 return gensizeerr(contextptr); 2443 p=cfs[1]; q=cfs[2]; r=cfs[0]; 2444 } else if (eq.type==_SYMB) { 2445 gen dy=identificateur(" dy"),d2y=identificateur(" d2y"),yx=symb_of(y,x); 2446 gen diffy=symbolic(at_derive,y),diff2y=symbolic(at_derive,diffy); 2447 eq=_subst(makesequence(eq,makevecteur(_derive(makesequence(yx,x),contextptr),_derive(makesequence(yx,x,2),contextptr)), 2448 makevecteur(dy,d2y)),contextptr); 2449 eq=_subst(makesequence(_subst(makesequence(_subst(makesequence(eq,diff2y,d2y),contextptr),diffy,dy),contextptr),yx,y),contextptr); 2450 if (eq.is_symb_of_sommet(at_equal)) 2451 eq=equal2diff(eq); 2452 eq=expand(eq,contextptr); 2453 vecteur terms=eq.is_symb_of_sommet(at_plus) && eq._SYMBptr->feuille.type==_VECT? 2454 *eq._SYMBptr->feuille._VECTptr:vecteur(1,eq); 2455 gen tmp; 2456 vecteur yvars=makevecteur(y,dy,d2y); 2457 for (const_iterateur it=terms.begin();it!=terms.end();++it) { 2458 if (is_constant_wrt_vars(tmp=_ratnormal(*it/y,contextptr),yvars,contextptr)) 2459 q+=tmp; 2460 else if (is_constant_wrt_vars(tmp=_ratnormal(*it/dy,contextptr),yvars,contextptr)) 2461 p+=tmp; 2462 else if (is_constant_wrt_vars(tmp=_ratnormal(*it/d2y,contextptr),yvars,contextptr)) 2463 r+=tmp; 2464 else return gensizeerr(contextptr); 2465 } 2466 } else return gensizeerr(contextptr); 2467 if (is_zero(r)) // not a second order ODE 2468 return gensizeerr(contextptr); 2469 p=_ratnormal(p/r,contextptr); 2470 q=_ratnormal(q/r,contextptr); 2471 if (rlvarx(p,x).size()+rlvarx(q,x).size()>2) // p or q is not rational in x 2472 return gensizeerr(contextptr); 2473 /* solve the equation y''+p(x)*y'+q(x)*y=0, transform it first to z''=r(x)*z */ 2474 r=sq(p)/4+_derive(makesequence(p,x),contextptr)/2-q; 2475 return kovacicsols(r,x,p,contextptr); 2476 } 2477 static const char _kovacicsols_s []="kovacicsols"; 2478 static define_unary_function_eval_quoted (__kovacicsols,&_kovacicsols,_kovacicsols_s); 2479 define_unary_function_ptr5(at_kovacicsols,alias_at_kovacicsols,&__kovacicsols,_QUOTE_ARGUMENTS,true) 2480 2481 /* 2482 * Return true iff the expression 'e' is constant with respect to 2483 * variables in 'vars'. 2484 */ is_constant_wrt_vars(const gen & e,const vecteur & vars,GIAC_CONTEXT)2485 bool is_constant_wrt_vars(const gen &e,const vecteur &vars,GIAC_CONTEXT) { 2486 for (const_iterateur it=vars.begin();it!=vars.end();++it) { 2487 if (!is_constant_wrt(e,*it,contextptr)) 2488 return false; 2489 } 2490 return true; 2491 } 2492 idnteval(const gen & g,GIAC_CONTEXT)2493 gen idnteval(const gen &g,GIAC_CONTEXT) { 2494 if (g.type==_IDNT) 2495 return _eval(g,contextptr); 2496 if (g.type==_SYMB) { 2497 gen &feu=g._SYMBptr->feuille; 2498 if (feu.type==_VECT) { 2499 vecteur v; 2500 for (const_iterateur it=feu._VECTptr->begin();it!=feu._VECTptr->end();++it) { 2501 v.push_back(idnteval(*it,contextptr)); 2502 } 2503 return symbolic(g._SYMBptr->sommet,change_subtype(v,feu.subtype)); 2504 } 2505 return symbolic(g._SYMBptr->sommet,idnteval(feu,contextptr)); 2506 } 2507 return g; 2508 } 2509 2510 #endif 2511 2512 #ifndef NO_NAMESPACE_GIAC 2513 } // namespace giac 2514 #endif // ndef NO_NAMESPACE_GIAC 2515