kill(all)$

maxima:special:peturbation (=second order,nonlinear bvp),version 1.0,part1
peter.vlasschaert@gmail.com,05/03/2020

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

load(lrats)$

packages : ODE = Ordinary Differential Equations

load(contrib_ode)$

Mathematical Model : 
===============

BVP = Boundary Value Problem 
----------------------------------------

p1:'diff(y,x,2)+y= cos(x)/(2+2*y^2);

p1a1:contrib_ode(p1,y,x);

p1a2:ode2(p1,y,x);

BC : Boundary Condition
------------------------------

1e) y(0) =2, x=0

p1b:y1 =2;

2e) y(π/2) = 1,x = π/2

p1c:y2=1;

first step : introduction ε
==================

p2:lhs(p1) =ε*rhs(p1);

find : zero-order problem ( ε=0)

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

 => find solution (general solution)

p2a1:ode2(p2a,y,x);

 => use BC 

p2aa:subst(x=0,rhs(p2a1))=rhs(p1b);
p2ab:subst(x=π/2,rhs(p2a1))=rhs(p1c);

 => solution (ε =0)

p2c:subst([p2aa,p2ab,y=y[0](x)],p2a1);

second step :perturbation series
========================

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

general : m → ∞

set m =2, second order approximation

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

first order :approximation
==================

p3b1:p3b,simpsum;

p3b2:subst([y=rhs(p3b1),ε^2=0],p2);

insert : zero -order approximation

p3b3:subst([p2c],p3b2);

p3b4:expand(ev(p3b3,nouns));

rem : y[0](x) ^2
============

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

expand((rhs(p2c))^2);

trigreduce(%);

expand(2*%)+2;

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

p4:denom(rhs(p3b2));

p4a:p4*lhs(p3b3) = num(rhs(p3b3));

p4b:expand(p4a);

p4c:ev(p4b,nouns);

p4d:subst([p2c],p4c);

p4e:expand(p4d);

p4f:coeff(p4e,ε,1);

p4g:trigsimp(p4f);

find coef : y[1](x)
----------------------
from lhs(p4g)

p4h:lhs(p4g);

p4h1:ratcoeff(p4h,y[1](x));

p4h2:trigreduce(p4h1);

p4h3:expand(p4h2);

final result:diff eq
=============

p5:part(p4h,1)+y[1](x)*p4h1 =  rhs(p4g);

see 'above'
--------------

p5a:lratsubst([p4h1=p4h3],p5);

final result : 'first order :approximation'

p5b:ratsimp(p5a/p4h3);

how to solve :ode2, p5b;

1e) y[1](x) = y
-------------------

p5b1:subst([y[1](x) =y],p5b);

y1=y[1](x)
example :variation parameter f(x) =perturb example,version 1.0,part2
peter.vlasschaert@gmail.com,03/03/2020

=>  y[n](0) = y[n](π/2) = 0 for n≥1

first order : solution
==============

result: y(x)=y[0](x) +y[1](x)
==========================
=>  y[0](0) = 2 ,y[0](π/2)
=>  y[n](0) = y[n](π/2) = 0 for n≥1
==========================