微分方程在建模中的应用（建立差分模型）

微分方程在建模中的应用（建立差分模型）
这里我们将k的值改变了下，这样更加符合算法的时间复杂度关系；
在温度扩散和浓度扩散模型中，微分方程形式为y'+g(x)*y=g(x)*f(x),目前我们假定g(x),f(x)均为常量，并且有g(x)=k,f(x)=m,下面解答这个方程并验证；
y'+k*y=k*m  (k,m为常量)
GO
采用标准解法有：
{ C*e^(kx) *  y}'=C*e^(kx)*k*m
GO
y=A+B*e^(-kx)
这里我们把k改为-1,
假设初始值y(0)=1,y(1)=10,k=-1
1=A+B
10=A+B*e^(1)
Go
A=1-B
B=9/(e-1)
注意将y=A+B*e^(-kx)回代入原方程中有ka=km,所以m=a


前面在用数值解法计算的时候采用的直接计算法，下面我们建立差分模型后直接计算它的值；
{y(n)-y(n-1)}/0.1+k*y(n-1)=k*m
Go
{y(n)-y(n-1)}/0.1-y(n-1)=-m
Go
y(n)-1.1*y(n-1)=-0.1*m
Go
根据差分方程的性质有：
y(n)=G*1.1^(10*n)+H  (其中乘以10的原因是为了它的每一步是微分方程对应一步的10倍)
这样有：
1=G+H
10=G*1.1^10+H
GO
G=9/(1.1^10-1)
H=1-G






下面写程序来证明：
(setq  k  -1)
(setq  b  (/ 9
             (- e  1)))
(setq  a  (- 1 b))
(setq  m  a)
(setq  g  (/ 9
             (- (pow 1.1 10)
                1)))
(setq  h  (- 1 g))


(defun pow (num count)
(if (or (> count 1) (eq  count  1) )
      (* num 
         (pow num 
              (- count 1) ) )
      1))


(defun slayer ( count)
(if (or (> count 1) (eq  count  1) )
      (* count 
         (slayer  
              (- count 1) ) )
      1))


(defun slayerex (num count)
(if (or (> count 1) (eq  count  1) )
      (* num 
         (slayerex 
              (-  num  1) 
              (- count 1) ) )
      1))




(defun  expr (x1 y1 )
(+  y1
    (*  1
        (-  (* k m)
            (* k y1)))))


(defun  calc (n x)
(if  (eq  n  0)
       1.0
    (+ (calc (1- n)
             x)
       (*  (pow  x  
                 n)
           (/    1
                 (slayer n))))))


(setq  e  (calc 10  1))
     


(defun  formula  (x)
(+   (*  b
         (pow e x))
     a))




(defun  formulaex  (x)
(+   (*  g  
         (pow  1.1
               (*  10  x)))
     h))




(defun  exprhelp (x1 )
(if (< (abs (-  x1  0))
       0.1)
        (formula  0)
    (expr  x1  
           (exprhelp  (-  x1  
                          1)))))                  


(defun  test (n)
(if (or (>  n 0) (eq  n  0))
  (progn 
       (print (formulaex n))
       (print  'compare)
       (print (formula n))      
       (test (- n 1)))
  (print 'over)))
[80]> (test  10)


77815.766
COMPARE
115365.836
29998.488
COMPARE
42438.04
11562.857
COMPARE
15609.402
4455.126
COMPARE
5739.699
1714.7881
COMPARE
2108.8384
658.26953
COMPARE
773.11945
250.93594
COMPARE
281.73593
93.8912
COMPARE
100.96606
33.34369
COMPARE
34.464542
10.0
COMPARE
10.0
1.0
COMPARE
1.0
OVER
OVER
注意因为这里是指数增长的，所以微分10步以内误差就比较大了，下面我们考虑继续减小差分的步长来缩短差距；新差分模型如下：
{y(n)-y(n-1)}/0.01+k*y(n-1)=k*m
Go
{y(n)-y(n-1)}/0.01-y(n-1)=-m
Go
y(n)-1.01*y(n-1)=-0.01*m
Go
根据差分方程的性质有：
y(n)=G*1.01^(100*n)+H  (其中乘以100的原因是为了它的每一步是微分方程对应一步的100倍)
这样有：
1=G+H
10=G*1.01^100+H
GO
G=9/(1.01^100-1)
H=1-G


(setq  g  (/ 9
             (- (pow 1.01 100)
                1)))
(setq  h  (- 1 g))


(defun  formulaex  (x)
(+   (*  g  
         (pow  1.01
               (*  100  x)))
     h))


[91]> (test  10)


110641.86
COMPARE
115365.836
40902.855
COMPARE
42438.04
15119.563
COMPARE
15609.402
5587.179
COMPARE
5739.699
2062.948
COMPARE
2108.8384
759.9984
COMPARE
773.11945
278.28305
COMPARE
281.73593
100.18735
COMPARE
100.96606
34.3433
COMPARE
34.464542
10.0
COMPARE
10.0
1.0
COMPARE
1.0
OVER
OVER
从这里可以看出来误差又小了点，如果我们继续这样的方法，显然就可以达到与微分解法几乎一样的结果；
在这里我们从理论上分析下原因：
y(n)=G*1.1^(10*n)+H这个差分结果我们可以做如下的变形：
Go
y(n)=G*(1+1/s)^(s*n)+H
Go
y(n)=G*{(1+1/s)^s}^n+H
Go
当s接近无穷大的时候(1+1/s)^s=e,所以有：
y(n)=G*e^n+H
这样就变为微分解法了；