kill(all)$

peter.vlasschaert@gmail.com, 16/12/2016

 part : intro

discrete : first-order nonhomogeous dynamical system
----------------------------------------------------
amount pollution => simplification to understand principles.

Lake Ontari:a(n)
Lake Eri   :b(n)

system of equations :difference equations

           a(n+1)=0.87*a(n)+0.38*b(n)    '1'
           b(n+1)=0.62*b(n)              '2'
 
solution : a(k) = c1*(0.87)^k+c2*(0.62)^k
           a(k) = general solution + particular solution
           rem : k= number of years > 0
           ?     a(k)
           ?     How to find c1,c2
           ?     How to get  eqn
pollutionlakemaximamodelsolvefinal.wxmx

eqn:4.52*(0.87)^k-1.52*(0.62)^k-0.3=0;

eqn: not 'polynomial' 

soln:allroots(eqn);

differentiation :   d(eqn)/dk

deqn  :'diff(eqn,k);
dkeqn : diff(eqn,k);

find zeros : eqn, ( for k > 0 ,year positiv number always)
=> close plot2d x

plot2d(part(eqn,1),[k,0,100])$

nu_f:ev(part(eqn,1),k=19.472);

number of years :Lake Ontario becomes : unpolluted

find_root(eqn,k,0,100);

algo:Newton-Raphson Method (example)
------------------------------------
func   :  equation 
x0     :  initial value
xf     :  final value
step_  :  steps
var    :  variable ,k
eps    :  error
   -----------------------------------------
   math : x(i+1) = x(i)-f(x(i))/(df(x(i)/dx)
   prog : xi     = xi - func/dfunc
   -----------------------------------------

one value    => (loop,go)

newton(func,var,x0,eps):=
block([xi,dfunc],
dfunc:diff(func,var),
xi:x0,
loop,
    xi:xi-subst(xi,var,func)/subst(xi,var,dfunc),
    if abs(subst(xi,var,func))<eps then return(xi),
go(loop)
    )$

func:lhs(eqn);

newton(func,k,19.0,10^(-4));

range values => (for,do)

newtonadv(func,var,x0,xf,step_,eps):=
block([xi,dfunc],
dfunc:diff(func,var),
for xi:x0 thru xf step step_ do 
  (    
    xi:xi-subst(xi,var,func)/subst(xi,var,dfunc),    
    if abs(subst(xi,var,func))<eps then return(xi)
  )   
 )$

rem: don't place ',' after return(xi),you get a error.

we want positive value,years,k>0

newtonadv(func,k,0,20,0.0001,10^(-4));

newtonadv(func,k,1,20,0.0001,10^(-4));

newtonadv(func,k,1.4,20,0.0001,10^(-4));