\( \DeclareMathOperator{\abs}{abs} \newcommand{\ensuremath}[1]{\mbox{$#1$}} \)
(%i1) kill(all)$

[email protected], 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

(%i1) eqn:4.52*(0.87)^k-1.52*(0.62)^k-0.3=0;
\[\tag{eqn}\label{eqn}4.52{{0.87}^{k}}-1.52{{0.62}^{k}}-0.3=0\]

eqn: not 'polynomial'

(%i2) soln:allroots(eqn);
\[\mbox{}\\\mbox{allroots: expected a polynomial; found }4.52{{0.87}^{k}}-1.52{{0.62}^{k}}-0.3\mbox{ -- an error. To debug this try: debugmode(true);}\]

differentiation :   d(eqn)/dk

(%i4) deqn  :'diff(eqn,k);
dkeqn : diff(eqn,k);
\[\tag{deqn}\label{deqn}\frac{d}{dk}\left( 4.52{{0.87}^{k}}-1.52{{0.62}^{k}}-0.3=0\right) \] \[\tag{dkeqn}\label{dkeqn}0.7266144174333597{{0.62}^{k}}-0.6294645443474546{{0.87}^{k}}=0\]

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

(%i5) plot2d(part(eqn,1),[k,0,100])$
(%i6) nu_f:ev(part(eqn,1),k=19.472);
\[\tag{nu\_ f}\label{nu\_ f}9.442236950424343{{10}^{-5}}\]

number of years :Lake Ontario becomes : unpolluted

(%i7) find_root(eqn,k,0,100);
\[\tag{\%{}o7}\label{o7} 19.47426223504164\]

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)

(%i8) 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)
   )$
(%i9) func:lhs(eqn);
\[\tag{func}\label{func}4.52{{0.87}^{k}}-1.52{{0.62}^{k}}-0.3\]
(%i10) newton(func,k,19.0,10^(-4));
\[\tag{\%{}o10}\label{o10} 19.47424609727112\]

range values ⇒ (for,do)

(%i11) 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

(%i12) newtonadv(func,k,0,20,0.0001,10^(-4));
\[\tag{\%{}o12}\label{o12} -3.086601656930104\]
(%i13) newtonadv(func,k,1,20,0.0001,10^(-4));
\[\tag{\%{}o13}\label{o13} \mathit{done}\]
(%i14) newtonadv(func,k,1.4,20,0.0001,10^(-4));
\[\tag{\%{}o14}\label{o14} 19.47419607178667\]
Created with wxMaxima.