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

update1: [email protected],19/12/2016

general equation (second-order linear equations)

discrete : second-order linear equation( ⇒ alternate approach)
-------------------------------------------------------------
p1*a(n+2)+p2*a(n+1)+p3*a(n)=0


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)\]

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)\]
(%i11) eqn_1:part(solve(eqn1,b(n)),1);
\[\mbox{}\\\mbox{rat: replaced -0.87 by -87/100 = -0.87}\mbox{}\\\mbox{rat: replaced -0.38 by -19/50 = -0.38}\] \[\tag{eqn\_ 1}\label{eqn\_ 1}\operatorname{b}(n)=\frac{100\operatorname{a}\left( n+1\right) -87\operatorname{a}(n)}{38}\]
(%i12) eqn_2:ev(eqn_1,n=n+1);
\[\tag{eqn\_ 2}\label{eqn\_ 2}\operatorname{b}\left( n+1\right) =\frac{100\operatorname{a}\left( n+2\right) -87\operatorname{a}\left( n+1\right) }{38}\]
(%i13) eqn_3:ev(eqn2,eqn_1,eqn_2);
\[\tag{eqn\_ 3}\label{eqn\_ 3}\frac{100\operatorname{a}\left( n+2\right) -87\operatorname{a}\left( n+1\right) }{38}=0.01631578947368421\left( 100\operatorname{a}\left( n+1\right) -87\operatorname{a}(n)\right) \]

p1*a(n+2)+p2*a(n+1)+p3*a(n)=0

(%i14) eqn_4:lhs(eqn_3)-rhs(eqn_3)=0;
\[\tag{eqn\_ 4}\label{eqn\_ 4}\frac{100\operatorname{a}\left( n+2\right) -87\operatorname{a}\left( n+1\right) }{38}-0.01631578947368421\left( 100\operatorname{a}\left( n+1\right) -87\operatorname{a}(n)\right) =0\]
(%i15) (0.62/0.38)/100;
\[\tag{\%{}o15}\label{o15} 0.01631578947368421\]

a(n+2) = x^2
a(n+1) = x
a(n)   = 1
p1*x^2+p2*x+p3 = 0, solve ?

(%i16) eqn_5:ev(eqn_4,a(n)=1,a(n+1)=x,a(n+2)=x^2);
\[\tag{eqn\_ 5}\label{eqn\_ 5}\frac{100{{x}^{2}}-87x}{38}-0.01631578947368421\left( 100x-87\right) =0\]
(%i17) eqn_6:ratsimp((100*x^2-87*x)/38-0.01631578947368421*(100*x-87)=0);
\[\mbox{}\\\mbox{rat: replaced -0.01631578947368421 by -31/1900 = -0.01631578947368421}\] \[\tag{eqn\_ 6}\label{eqn\_ 6}\frac{5000{{x}^{2}}-7450x+2697}{1900}=0\]
(%i18) eqn_7:float(eqn_6);
\[\tag{eqn\_ 7}\label{eqn\_ 7}5.263157894736842{{10}^{-4}}\,\left( 5000.0{{x}^{2}}-7450.0x+2697.0\right) =0.0\]
(%i19) eqn_8:part(eqn_7,1,2);
\[\tag{eqn\_ 8}\label{eqn\_ 8}5000.0{{x}^{2}}-7450.0x+2697.0\]
(%i22) p1:coeff(eqn_8,x^2);
p2:coeff(eqn_8,x);
p3:ev(eqn_8,x=0);
\[\tag{p1}\label{p1}5000.0\] \[\tag{p2}\label{p2}-7450.0\] \[\tag{p3}\label{p3}2697.0\]
(%i23) eqs:solve(eqn_5,x);
\[\mbox{}\\\mbox{rat: replaced -0.01631578947368421 by -31/1900 = -0.01631578947368421}\] \[\tag{eqs}\label{eqs}[x=\frac{87}{100},x=\frac{31}{50}]\]
(%i24) 31.0/50.0;
\[\tag{\%{}o24}\label{o24} 0.62\]
(%i26) pc1:float(part(eqs,1,2));
pc2:float(part(eqs,2,2));
\[\tag{pc1}\label{pc1}0.87\] \[\tag{pc2}\label{pc2}0.62\]
(%i27) eq_model:a(k)=c1*(pc1)^k+c2*(pc2)^k;
\[\tag{eq\_ model}\label{eq\_ model}\operatorname{a}(k)=\mathit{c1}\,{{0.87}^{k}}+\mathit{c2}\,{{0.62}^{k}}\]

given : a(0),a(1) to find c1,c2

We want to find value for a(1),a(0) = fixed value = 3
a(0) = 3  ( Why a(0) = 3*b(0),b(0)=1,three times as large)

(%i29) eq_model1:ev(eq_model,k=0,a(k)=3);
eq_model2:ev(eq_model,k=1);
\[\tag{eq\_ model1}\label{eq\_ model1}3=1.0\mathit{c2}+1.0\mathit{c1}\] \[\tag{eq\_ model2}\label{eq\_ model2}\operatorname{a}(1)=0.62\mathit{c2}+0.87\mathit{c1}\]

solve the system : a(0),a(1)

(%i30) cf:solve([eq_model1,eq_model2],[c1,c2]);
\[\mbox{}\\\mbox{rat: replaced -1.0 by -1/1 = -1.0}\mbox{}\\\mbox{rat: replaced -1.0 by -1/1 = -1.0}\mbox{}\\\mbox{rat: replaced -0.87 by -87/100 = -0.87}\mbox{}\\\mbox{rat: replaced -0.62 by -31/50 = -0.62}\] \[\tag{cf}\label{cf}[[\mathit{c1}=\frac{100\operatorname{a}(1)-186}{25},\mathit{c2}=-\frac{100\operatorname{a}(1)-261}{25}]]\]
(%i32) cf1:part(cf,1,1);
cf2:part(cf,1,2);
\[\tag{cf1}\label{cf1}\mathit{c1}=\frac{100\operatorname{a}(1)-186}{25}\] \[\tag{cf2}\label{cf2}\mathit{c2}=-\frac{100\operatorname{a}(1)-261}{25}\]
(%i33) model_lake:ev(eq_model,cf1,cf2);
\[\tag{model\_ lake}\label{model\_ lake}\operatorname{a}(k)=\frac{\left( 100\operatorname{a}(1)-186\right) \,{{0.87}^{k}}}{25}-\frac{\left( 100\operatorname{a}(1)-261\right) \,{{0.62}^{k}}}{25}\]
(%i34) eqn_check :aa(k)=4.52*(0.87)^k-1.52*(0.52)^k;
\[\tag{eqn\_ check}\label{eqn\_ check}\operatorname{aa}(k)=4.52{{0.87}^{k}}-1.52{{0.52}^{k}}\]
(%i35) eqn_check_inp:ev(eqn_check,k=1);
\[\tag{eqn\_ check\_ inp}\label{eqn\_ check\_ inp}\operatorname{aa}(1)=3.142\]
(%i36) a_1:part(eqn_check_inp,2);
\[\tag{a\_ 1}\label{a\_ 1}3.142\]
(%i38) cf2a1:ev(cf2,a(1)=a_1);
cf1a1:ev(cf1,a(1)=a_1);
\[\tag{cf2a1}\label{cf2a1}\mathit{c2}=-2.127999999999997\] \[\tag{cf1a1}\label{cf1a1}\mathit{c1}=5.127999999999997\]
(%i39) eq_model_a1:ev(eq_model,cf2a1,cf1a1);
\[\tag{eq\_ model\_ a1}\label{eq\_ model\_ a1}\operatorname{a}(k)=5.127999999999997{{0.87}^{k}}-2.127999999999997{{0.62}^{k}}\]

check : a(0),a(1)

(%i40) eq_model_a_1:ev(eq_model_a1,k=1);
\[\tag{eq\_ model\_ a\_ 1}\label{eq\_ model\_ a\_ 1}\operatorname{a}(1)=3.141999999999999\]
(%i41) eq_model_a_0:ev(eq_model_a1,k=0);
\[\tag{eq\_ model\_ a\_ 0}\label{eq\_ model\_ a\_ 0}\operatorname{a}(0)=3.0\]
Created with wxMaxima.