(%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; |
eqn: not 'polynomial'
(%i2) | soln:allroots(eqn); |
differentiation : d(eqn)/dk
(%i4) |
deqn :'diff(eqn,k); dkeqn : diff(eqn,k); |
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); |
number of years :Lake Ontario becomes : unpolluted
(%i7) | 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)
(%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); |
(%i10) | newton(func,k,19.0,10^(-4)); |
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)); |
(%i13) | newtonadv(func,k,1,20,0.0001,10^(-4)); |
(%i14) | newtonadv(func,k,1.4,20,0.0001,10^(-4)); |