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

[email protected],18/02/2016

part : how to get the model

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

Lake Ontari:a(n)
Lake Eri   :b(n)
Lake Ontario connected with Lake Eri

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'

(%i2) eqn1:a(n+1)=a1*a(n)+b1*b(n)$
eqn2:b(n+1)=a2*a(n)+b2*b(n)$

b(n):each year: Lake Eri 38%(0.38) 'replaced' by rain

(%i4) a2:0;
b2:1-0.38;
\[\tag{a2}\label{a2}0\] \[\tag{b2}\label{b2}0.62\]
(%i6) eqn1:a(n+1)=a1*a(n)+b1*b(n);
eqn2:b(n+1)=a2*a(n)+b2*b(n);
\[\tag{eqn1}\label{eqn1}\operatorname{a}\left( n+1\right) =\mathit{b1}\,\operatorname{b}(n)+\mathit{a1}\,\operatorname{a}(n)\] \[\tag{eqn2}\label{eqn2}\operatorname{b}\left( n+1\right) =0.62\operatorname{b}(n)\]

a(n);each year :Lake Ontario
38%(0.38) Lake Eri flows into Lake Ontario
13%(.13 ) Lake Ontarion 'replaced' from water from Lake Eri

rem:replaced : 1-..

(%i8) a1:1-0.13;
b1:0.38;
\[\tag{a1}\label{a1}0.87\] \[\tag{b1}\label{b1}0.38\]
(%i10) eqn1:a(n+1)=a1*a(n)+b1*b(n);
eqn2:b(n+1)=a2*a(n)+b2*b(n);
\[\tag{eqn1}\label{eqn1}\operatorname{a}\left( n+1\right) =0.38\operatorname{b}(n)+0.87\operatorname{a}(n)\] \[\tag{eqn2}\label{eqn2}\operatorname{b}\left( n+1\right) =0.62\operatorname{b}(n)\]

solve: problem from the connected lakes(Eri,Ontario)

first: not connected  :eqn2 : n=0,1,2,3,...
b(0)=1 unit of pollution

(%i13) v0:b(0)=1;
v1:b(1)=0.62*b(0);
v2:b(2)=0.62*b(1);
\[\tag{v0}\label{v0}\operatorname{b}(0)=1\] \[\tag{v1}\label{v1}\operatorname{b}(1)=0.62\operatorname{b}(0)\] \[\tag{v2}\label{v2}\operatorname{b}(2)=0.62\operatorname{b}(1)\]
(%i14) w1:ev(v2,v1);
\[\tag{w1}\label{w1}\operatorname{b}(2)=0.3844\operatorname{b}(0)\]
(%i15) d1:part(coeff(w1,b(0)),2);
\[\tag{d1}\label{d1}0.3844\]
(%i17) d2:d1-(b2)^kk=0;
find_root(d2,kk,0,100);
\[\tag{d2}\label{d2}0.3844-{{0.62}^{\mathit{kk}}}=0\] \[\tag{\%{}o17}\label{o17} 2.0\]

b(k)=(0.62)^k*b(0)=(0.62)^k*1

(%i19) eqn11:ev(eqn1,b(n)=(0.62)^n*b(0));
eqn111:ev(eqn11,b(0)=1);
\[\tag{eqn11}\label{eqn11}\operatorname{a}\left( n+1\right) =0.87\operatorname{a}(n)+0.38\operatorname{b}(0)\,{{0.62}^{n}}\] \[\tag{eqn111}\label{eqn111}\operatorname{a}\left( n+1\right) =0.87\operatorname{a}(n)+0.38{{0.62}^{n}}\]

theory: Nonhomogeneneous dynamical systems
       a(n+1)=r*a(n)+b*s^n
suppose  : r<>s
solution : a(k)=c1*(r)^k+c2*(s)^k
          a(k)=c*(r)^k+a*s^k
          a(k)= general solution + particular solution
          ? find a
rem: r=s
    a(k)=c1*(r)^k+particular solution
    a(k)=c1*(r)^k+c2*k*(r)^k,c2=b/r
    a(k)=c1*r^k+(b/r*k)*r^k

(%i22) eqt : a(n+1)=r*a(n)+b*s^n;
eqak: a(k)=c*(r)^k+a*s^k;
eqan: ev(eqak,k=n);
\[\tag{eqt}\label{eqt}\operatorname{a}\left( n+1\right) =b\,{{s}^{n}}+\operatorname{a}(n)r\] \[\tag{eqak}\label{eqak}\operatorname{a}(k)=a\,{{s}^{k}}+c\,{{r}^{k}}\] \[\tag{eqan}\label{eqan}\operatorname{a}(n)=a\,{{s}^{n}}+c\,{{r}^{n}}\]
(%i23) eqanplus1:ev(eqan,n=n+1);
\[\tag{eqanplus1}\label{eqanplus1}\operatorname{a}\left( n+1\right) =a\,{{s}^{n+1}}+c\,{{r}^{n+1}}\]
(%i24) eqinput:ev(eqt,eqan,eqanplus1);
\[\tag{eqinput}\label{eqinput}a\,{{s}^{n+1}}+c\,{{r}^{n+1}}=r\,\left( a\,{{s}^{n}}+c\,{{r}^{n}}\right) +b\,{{s}^{n}}\]
(%i25) eqval1:ratsimp(a*s^(n+1)+c*r^(n+1)=r*(a*s^n+c*r^n)+b*s^n);
\[\tag{eqval1}\label{eqval1}a\,{{s}^{n+1}}+c\,{{r}^{n+1}}=\left( ar+b\right) \,{{s}^{n}}+c\,{{r}^{n+1}}\]
(%i26) radcan(a*s^(n+1)+c*r^(n+1)=(a*r+b)*s^n+c*r^(n+1));
\[\tag{\%{}o26}\label{o26} a\,{{s}^{n+1}}+c\,{{r}^{n+1}}=\left( ar+b\right) \,{{s}^{n}}+c\,{{r}^{n+1}}\]
(%i30) d_part1:coeff(part(eqval1,1),s^(n+1))*s;
d_part2:coeff(part(eqval1,2),s^(n));
eq_particular: d_part1=d_part2;
aa:part(solve(eq_particular,a),1);
\[\tag{d\_ part1}\label{d\_ part1}as\] \[\tag{d\_ part2}\label{d\_ part2}ar+b\] \[\tag{eq\_ particular}\label{eq\_ particular}as=ar+b\] \[\tag{aa}\label{aa}a=\frac{b}{s-r}\]
(%i31) gen_eq1:ev(part(eqak,2),aa);
\[\tag{gen\_ eq1}\label{gen\_ eq1}\frac{b\,{{s}^{k}}}{s-r}+c\,{{r}^{k}}\]
(%i32) gen_eq:a(k)=gen_eq1;
\[\tag{gen\_ eq}\label{gen\_ eq}\operatorname{a}(k)=\frac{b\,{{s}^{k}}}{s-r}+c\,{{r}^{k}}\]

example: from Lake Ontario and Lake Eri
        r=0.87
        s=0.62
        b=0.38
size: Lake Ontario three times the size of Lake Eri
result :three times pollution but the same initial concentration
       a(0)=a0 ,a0=3*b0=3.b(0)

(%i33) gen_eq_example:ev(gen_eq,r=0.87,s=0.62,b=0.38);
\[\tag{gen\_ eq\_ example}\label{gen\_ eq\_ example}\operatorname{a}(k)=c\,{{0.87}^{k}}-1.52{{0.62}^{k}}\]

because c undetermined variable,anything we like.
check initial conditions.

(%i35) gen_eq_ex1:ev(gen_eq_example,c=c-part(aa,2));
gen_eq_ex2:ev(gen_eq_ex1,r=0.87,s=0.62,b=0.38);
\[\tag{gen\_ eq\_ ex1}\label{gen\_ eq\_ ex1}\operatorname{a}(k)={{0.87}^{k}}\,\left( c-\frac{b}{s-r}\right) -1.52{{0.62}^{k}}\] \[\tag{gen\_ eq\_ ex2}\label{gen\_ eq\_ ex2}\operatorname{a}(k)=\left( c+1.52\right) \,{{0.87}^{k}}-1.52{{0.62}^{k}}\]
(%i36) gen_eq_init:ev(gen_eq_ex2,k=0);
\[\tag{gen\_ eq\_ init}\label{gen\_ eq\_ init}\operatorname{a}(0)=1.0\left( c+1.52\right) -1.52\]
(%i37) gen_eq_initc:ev(gen_eq_init,a(0)=3);
\[\tag{gen\_ eq\_ initc}\label{gen\_ eq\_ initc}3=1.0\left( c+1.52\right) -1.52\]
(%i38) c_init:part(solve(gen_eq_initc,c),1);
\[\mbox{}\\\mbox{rat: replaced 4.52 by 113/25 = 4.52}\mbox{}\\\mbox{rat: replaced -1.0 by -1/1 = -1.0}\mbox{}\\\mbox{rat: replaced 1.52 by 38/25 = 1.52}\] \[\tag{c\_ init}\label{c\_ init}c=3\]
(%i39) gen_eq_ex3:ev(gen_eq_ex2,c_init);
\[\tag{gen\_ eq\_ ex3}\label{gen\_ eq\_ ex3}\operatorname{a}(k)=4.52{{0.87}^{k}}-1.52{{0.62}^{k}}\]

How long will take :Lake Eri , 10%(0.1) from present value b(0)=1
Because Lake Ontario '3' times size of Lake Eri
                              30%(0.3)
(Lake Eri :ref)

first :unpolluted,Lake Eri ? How many years

(%i40) b_unpol1:ev(b(k)=(0.62)^k*b(0),b(0)=1);
\[\tag{b\_ unpol1}\label{b\_ unpol1}\operatorname{b}(k)={{0.62}^{k}}\]
(%i41) b_unpol2:ev(b_unpol1,b(k)=0.1);
\[\tag{b\_ unpol2}\label{b\_ unpol2}0.1={{0.62}^{k}}\]

find ,between 0 years and 100 years

(%i42) k_year_unpol_Eri:find_root(b_unpol2,k,0,100);
\[\tag{k\_ year\_ unpol\_ Eri}\label{k\_ year\_ unpol\_ Eri}4.816762862638821\]

second:unpolluted,Lake Ontario ? How many years
b(k) = 0.1 → 3 times → a(k) = .3

(%i43) a_unpol1:ev(gen_eq_ex3,a(k)=.3);
\[\tag{a\_ unpol1}\label{a\_ unpol1}0.3=4.52{{0.87}^{k}}-1.52{{0.62}^{k}}\]

find,between 0 years and 100 years

(%i44) k_year_unpol_Ontario:find_root(a_unpol1,k,0,100);
\[\tag{k\_ year\_ unpol\_ Ontario}\label{k\_ year\_ unpol\_ Ontario}19.47426223504164\]
Created with wxMaxima.