kill(all)$

maxima: intro to solve with python , ode(finite difference),version 1.0a
peter.vlasschaert@gmail.com,31/08/2022

load packages:

load(lrats)$ ratprint:false$

math :simple ode : u"+w^2*u = 0 (linear second order differential equation)
ic 1:u'(0) = r
ic 2:u(0)  =1
ic = (i)nitial (c)ondition

ode : (o)rdinary (d)ifferential (e)quation

→ exact solution : maxima

p1:'diff(u,t,2)+w^2*u =0;

assume(w>0)$

p1a:ode2(p1,u,t);

ic : find {%k1,%k2}

ic1 → u'(t) = u'(0) =0

p1b:'diff(u,t) = diff(rhs(p1a),t);

p1b1:subst([t=0],factor(rhs(p1b))) = 0; 

assume : w>0

p1b1a:[part(lhs(p1b1),1) =rhs(p1b1)];

result : p1a

p1a1: u =subst(p1b1a,rhs(p1a));

ic2 → u(t) = u(0) =1 ,use : p1a1 → u = 1,t=0

p1a2:[reverse(subst([u=1,t=0],p1a1))];

final result : ode 

p1a3:subst(p1a2,p1a1);

check :in ode with solution ⇒ p1a3
→ Answer = true or false

p1a4:subst([p1a3],p1);

p1a41:is(ev(p1a4,nouns));

→ finite difference ⇒ (use central difference = u'') ⇒ f(U(.))
rem : . = {n+1,n,n-1},see : below

1e step:u''= 'diff(u,t,2)

p2a:part(lhs(p1),2) =(U(n+1)-2*U(n)+U(n-1))/(Δt^2);

2e step :u' = 'diff(u,t)

p2b:'diff(u,t) = (U(n+1)-U(n-1))/(2*Δt);

3e step :u = U(n)

p2c: u = U(n);

4e step: final : difference equation : from"ode"

p3:subst([p2a,p2b,p2c],p1);

5e step : ic : n=0, (t=0,ode)

p3a:subst([n=0],p3);

6e step : solve U(1)

p3b:solve(p3a,U(1))[1];

7e step : How to find  : U(-1) , -1 means -Δt

solve : find from p2b ⇒with n=0 and ic 1

p3c:subst([n=0],rhs(p2b)) =r;/*simple r=0*/;

p3c1:solve(p3c,U(-1));

8e step : before final  result when :subst p3b ,use p3c1

p4:subst(p3c1,p3b);

9e step : before final result (use ic2 , u(0) =1), use p2c

p4a:reverse(subst([n=0,u=1],p2c));

10e step: final result :"difference equation"

p5:subst([p4a],p4);

when r = 0 for u'(0) =r,special case, general : float

p5d:subst([r=0],p5);

when r ≠ 0 for u'(0) =r,r : float

p5a:ratsimp(solve(p5,U(1)));

special case : r=0 for u`(0) =0

p5d1:subst([r=0],p5a);

11e step : find ⇒ U(n+1),p3, recursion difference equation.

p5b:solve(p3,U(n+1))[1];

rem : express this in python , N steps

example : n=1
----------------

p5c:subst([n=1],p5b);

example : n=N, result: "N+1 points in the domain"
-----------------------------------------------------------

p5c1:subst([n=N],p5b);

----------------------------------------------------------