kill(all)$

Maxima 'special topics' : Introduction to perturbations   version 1.0
peter.vlasschaert@gmail.com,02/03/2020

load packages
===========

load(lrats);

rem : one variable problems

def :Perturbation theory ( = collection of methods) to get approximate solution
       to problems involving (small parameter = ε).

def :singular perturbation 
==================
Whose solution ε = 0, different in character from "neighbourhood",
solution obtainded in the (limit , ε  → 0).

Regular perturbation
===============

a) roots of polynomial
================

p1:x^5-16*x+1 = 0;

a1) find roots 'p1'
---------------------

p1a:allroots(p1);

alll x roots.'list'
----------------

p1a1:makelist(rhs(p1a[i]),i,1,5);

result : five roots expected = 3 real roots + 2 complex roots (*complex conjugate =z*).
==============================================================
complex conjugate 
==============
z  = a+i*b  => z* = a-i*b

a2) how to do introduction  of ε non perturbation problem
---------------------------------------------------------------------
not perturbation problem = 'p1'

try : to replace constant term,'p1' with ε

p2:ev(lhs(p1),x=0)=ε;

p2a:lratsubst([p2],p1);

a3) what relation needed to connect ε with x
================================

finite
====

p3a:x(ε,m) = sum(a[n]*ε^n, n,0,m);

infinite
=====

p3b:lhs(p3a)=limit(rhs(p3a),m,inf);

p3b1:lhs(p3a) = subst([m=inf],rhs(p3a));

result : m → ∞

p3b2:rhs(p3b) =rhs(p3b1);

sol1 :how first term  :serie 'p2a' with ε =0

p4:subst([ε = 0],p2a);

p4a:gfactor(p4);

p4b:allroots(p4a);

 => zeroth-order 'perturbation theory'

p4c:makelist(rhs(p4b[i]),i,1,5);

p4d:float(p4c[3]);

p4da:-2;

rem : five roots ,we want approximate 3rd root

 => second-order 'perturbation theory' , 'p4d'

we use : 'finite' , m=2 into 'p3a'

p4d1:subst([m=2],p3a);

p4d1a:p4d1,simpsum;

p4d1b:subst([a[0]=p4da],p4d1a);

 => subst 'above' into p2a
 ε^k = 0 , k=3,......
why : only two unknows,a[1],a[2]

p4e:subst([x=rhs(p4d1b)],p2a);

p4e1:expand(p4e);

p4e1a:makelist(ε^k=0,k,3,12);

p4e1b:subst(p4e1a,p4e1);

factor from ε^2
===========

1e) equation

p4e1b1:coeff(p4e1b,ε,2);

factor from ε
===========

2e) equation

p4e1b2:coeff(p4e1b,ε,1);

solve : 1e) & 2e) equation

p4f:solve([p4e1b1,p4e1b2],[a[1],a[2]]);

p4f1:flatten(p4f);

first :approximation , root = -2
======================
use : p4d1b,with ε = 1

p5:subst(p4f1,p4d1b);

limit : ε → 0

p5:float(subst([ε=1],p5));

check : -2,'p5'
===========
orginal equation,'p1'

p6a:subst([x=-2],p1);

bad :solution ( = very bad)

p6b:subst([x=rhs(p5)],p1);

good : solution

four steps: perturbative analysis
========================

1e)  step :introduction ε, 'orginal equation,here constant polynomial equation '
2e)  introduction perturbation series ,see,'p3a' & 'p3b', (example : ε^k=0 for 3 to max(powerof ε))
3e)  answer:find coef : a[i] (example : m=2,a[1],a[2], ε→1), x(ε,m) =x(1,m)
4e) check with orginal equation ,without ε.

rem :a[0] = root in question ( a[0] =-2.0)