kill(all)$

maxima :point collocation  method,version 1.0,part 2
peter.vlasschaert@gmail.com,17/11/2020

see:website
maxima :point collocation  method,version 1.0,part 1

load packages , option.
==================

load(lrats)$ ratprint:false$

powerdisp:true$

part 2 : not simple polynoom
======================

G =G(x,y,dy/dx,d/dx(dy/dx),........)

rem :G =G(x,y,dy/dx)  ,part : 1

problem :
=======
 →  ODE :d(dy/dx)/dx -y +x=0, domain : 0≤x≤1,
 "begin :x=0 → y(0) =0  & end: x=1 → y(1) =0",domain : 0≤ x ≤ 1

p1:'diff(y,x,2)-y+x =0;

p11:p1-x;

1e) step :approximate solution 
------------------------------------------
a) 2e order nonlinear ODE (BVP  : Boundary Value Problem:initial & final)
b)  approximation second order polynoom : 'not simple polynoom'

p1a:y=sum(a[i]*sin(i*%pi*x),i,1,3);

p1a1:ev(p1a,numer);

2e) BC 
------------

→  BVP 'begin'  :  x= 0 → y=0
------------------------------------------

p1b:subst([x=0,y=0],p1a);

p1b1:[reverse(p1b)];

→  BVP 'end'  :  x = 1→ y=0
--------------------------------------------

p1b2:subst([x=1,y=0],p1a);

rem :example here are BCs inside : trial function ( p1a,y)

3e)  residual  = R = 0 
------------------------------
---------------------------------------------------------------------------------------------------
find:  a[i] , i=1,2,3
---------------------------------------------------------------------------------------------------
a) R=0   → x=0.25 , value inside interval : 0 ≤ x =0.25 ≤ 1
b) R=0   → x=0.50 , value inside interval : 0 ≤ x =0.50 ≤ 1
c) R=0   → x=0.75 , value inside interval : 0 ≤ x =0.75 ≤ 1
----------------------------------------------------------------------------------------------------
reason : a[i], i=1,2,3  => need three collocation points ' solve the problem'

p1c:R = lhs(p1);

p1ca:subst([p1a],p1c);

p1ca:ev(p1ca,nouns);

a) R=0   → x=0.25

p1c1b1:ev(subst([R=0,x=0.25],p1ca),nouns),numer;

b) R=0   → x=0.50

p1c1b2:ev(subst([R=0,x=0.50],p1ca),nouns),numer;

b) R=0   → x=0.75

p1c1b3:ev(subst([R=0,x=0.75],p1ca),nouns),numer;

4e)  find a[i], i=1,2,3
----------------------------

p1d:listofvars(p1c1b1);

p1d1:solve([p1c1b1,p1c1b2,p1c1b3],p1d)[1],numer;

5e) subst  'p1d' into p1a1
------------------------------------

p1d2:subst(p1d1,p1a1);

p1d2a:float(p1d2);

→  approximation  solution : ODE 

wxplot2d([rhs(p1d2)], [x,0,1])$

find :exact solution
---------------------------

→  use ode2

p1e:ode2(p1,y,x);

p1ea:listofvars(p1e);

→  find : %k1,%k2

p1e1:subst([x=0,y=0],p1e);

p1e1a:[reverse(p1e1)];

p1e1b:subst([x=1,y=0],p1e);

→  exact solution

p1e1c:flatten(solve([p1e1,p1e1b],p1ea));

p1e1ca:[p1e1c[3],p1e1c[4]];
p1e1cb:subst(p1e1ca,p1e);

wxplot2d([rhs(p1e1cb)], [x,0,1])$

see : final result on the graph below
===========================

wxplot2d([rhs(p1d2),rhs(p1e1cb)], [x,0,1])$

rem :approximation very good one see wxplot2d,not always know exact solution.
-------------------------------------------------------------------------------------------------------------
blue : approximation 
red   : exact function
--------------------------------------------------------------------------------------------------------------