kill(all)$

math:solve differential equations with infinite sums( intro),version 1.0
peter.vlasschert@gmail.com,17/04/2020

load ,packages
=============

load(simplify_sum)$

load(lrats)$

How to changevar , "Maxima"
=======================

'integrate(sin(3*x),x,0,4);

changevar (%, y=3*x, y, x);

sum (a[i]*x^(i-2), i, 0, inf);

changevar (%, i=n+2,n,i);

differential equations:  y'-y=0 (1)
======================
rem : dy/dx =y'(x), y=y(x)

1e) We assume form 
=================

y=a0+a1*x+a2*x^2+a3*x^3+......... (2)

Mathematical : representation 'Maxima'
=======================

p1:y=sum (a[i]*x^i, i, 0, inf);

1e) step : find  y'
-----------------------

p1a:'diff(y,x) =diff(rhs(p1),x);

2e) write differential equation (1)
---------------------------------------------

p2a:p1a-p1;

result : differential equations 
---------------------------------------

p2a1:lhs(p2a)=0;

p2a2:rhs(p2a)=0;

result : p2a1  → (1)
            p2a2  → approximation with sums

3e)  find coef x^n
-------------------------

p2a3:subst([i=n],p2a2);

p2a4:part(lhs(p2a2),1);

p2a5:changevar (%, i=n+1,n,i);

p2a6:part(lhs(p2a3),2);

p2a7:intosum(p2a6);

to do an 'intosum (expr)' before the 'sumcontract'. !!!!!

p2a8:sumcontract(p2a5+p2a7)=rhs(p2a3);

find factor before :x^n
------------------------------
math : factor = 0 (3)

p2a9:factor(p2a8);

lhs : 'p2a9'

p2a10:args(lhs(p2a9));

see:(3),above.

p2a11:ratcoef(p2a10,x^n)[1]=0;

find:a[n+1]
----------------
use:p2a11

p2a12:solve(p2a11,a[n+1])[1];

n : integer.
--------------

p3:makelist(ev(p2a12,nouns),n,0,7);

p3a:length(p3);

option :change
============
→ factorflag ?

default: option

factorflag;

change option

factorflag:true;

why this option :factor a number in primes

find :a[2]=f(a[0])
---------------------

even : 'n'
=======

p3b:factor(subst(p3[1],p3[2]));

find:a[3]=f(a[0])
---------------------

odd: 'n'
======

p3c:factor(subst(p3b,p3[3]));

from here : a[n]=a[0]/n!

How to find : a[8]=f(a[0])
-----------------------------------

p4:p3;

p4a:lhs(p4[1])=rhs(p4[1]);

p4b:lhs(p4[2])=factor(subst(%,rhs(p4[2])));

p4c:lhs(p4[3])=factor(subst(%,rhs(p4[3])));

p4d:lhs(p4[4])=factor(subst(%,rhs(p4[4])));

p4e:lhs(p4[5])=factor(subst(%,rhs(p4[5])));

p4f:lhs(p4[6])=factor(subst(%,rhs(p4[6])));

p4f:lhs(p4[7])=factor(subst(%,rhs(p4[7])));

p4f:lhs(p4[8])=factor(subst(%,rhs(p4[8])));

denom(rhs(%));

check for 8!
-----------------

float(8!);

use : possibble automatization

find solution :a[n] for 
================

p5:a[n]=a[0]/n!;

p5a:subst([i=n,p5],p1);

to find  ?  p5a
===========

p5b:lhs(p5a)=args(rhs(p5a));

p5c:rhs(p5b)[2],simplify_sum;

final solution
==========

p5d:lhs(p5b)=rhs(p5b)[1]*p5c;

check :solution (= back substitution)
============
use :p2a

a) subst
-----------

p6:subst([p5d],p2a1);

b) (ev)ulation, nouns
-----------------------------

p6a:ev(p6,diff);

c) boolean 
---------------
is, command 'Maxima'

p6b:is(p6a);

check solution :ode2
-----------------------------
use:p1a,p1

a) ode to solve

p7:lhs(p1a-p1)=0;

b) solution

p7a:ode2(p7,y,x);

c) find :%c
use:p7a,p5d

p7b:rhs(p7a)=rhs(p5d);

p7c:solve(p7b,%c);

finish:

p8:p4;

rhs(p8[1]);

p9:p2a12;

p9a:makelist(ev(lhs(p9),nouns),n,0,7);

p9b:makelist(ev(rhs(p9),nouns),n,0,7);

p9c:map("+",p9b,p9a);

p10c:p8;

ev(lratsubst([a(k-1)=a(k-2)/(k-1)],a(k)=a(k-1)/k),nouns)

m(k):=lratsubst([a(k-1)=a(k-2)/(k-1)],a(k)=a(k-1)/k);

p10d:create_list(m(j),j,[2,3,4,5,6,7,8]);

q:1;

p10e:subst([p10d[q]],p10d[q+2]);

p10f:makelist(subst([p10d[q]],p10d[q+2]),q,1,length(p10d)-2);

last step:see 'p4f'

p10g:subst([first(p10f)],last(p10f));

p10h:makelist(subst([p8[q]],p8[q+1]),q,1,length(p8)-1);

p10h1:subst([p8[1]],p10h);

p10h2:makelist(lratsubst([p10h1[q]],p10h1[q+1]),q,1,length(p10h1)-1,2);

p10h3:makelist(lratsubst([p10h2[q]],p10h2[q+1]),q,1,length(p10h2)-1);

p10h4:subst([first(p10h2)],last(p10h3));

? 5040

p10h5:float(part(rhs(p10h4),2));

p10h6:float(7!);

===========================================================================