伯努利方法的应用(也是自治微分方程)
假设微分方程模型为y'=y+1/y
方程两边分别除以1/y有:
y'/{1/y}=y/{1/y}+1
Go
y'/y^(-1)=1/y^(-2)+1
按照伯努利方法设置v=1/y^(-2),则v'=2*y'*/y^(-1),
GO
v'/2=v+1
Go
v'=2v+2
Go
v'-2v=2
左边部分按照标准解法有:
|u|=C*e^(-2x)
GO
可以得到v=-1/2+C*e^(2x)
这样可以得到|y|=sqrt v =sqrt {-1/2+C*e^(2x)} ,
现在假设初始条件是y(0)=3,那么C=9.5
¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥
下面写程序来证明:
(setq c 9.5)
(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))
y'=y+1/y
(defun expr (x1 y1 )
(+ y1
(* 0.01
(+ y1
(/ 1
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)
(sqrt
(+ -0.5
(* c
(pow e
(* 2
x))))))
(defun exprhelp (x1 )
(if (< (abs (- x1 0))
0.01)
(formula 0)
(expr x1
(exprhelp (- x1
0.01)))))
(defun test (n)
(if (> n 0)
(progn
(print (exprhelp n))
(print 'compare)
(print (formula n))
(test (- n 1)))
(print 'over)))
[25]> (test 20)
1.3757553E9
COMPARE
1.4953814E9
5.08632E8
COMPARE
5.5012006E8
1.8804696E8
COMPARE
2.0237784E8
6.9523024E7
COMPARE
7.445064E7
2.5703448E7
COMPARE
2.7388858E7
9502851.0
COMPARE
1.0075797E7
3513310.5
COMPARE
3706678.5
1298910.3
COMPARE
1363610.8
480221.44
COMPARE
501644.34
177543.2
COMPARE
184544.63
65639.7
COMPARE
67890.164
24267.734
COMPARE
24975.395
8972.051
COMPARE
9187.934
3317.067
COMPARE
3380.0518
1226.3562
COMPARE
1243.4513
453.3967
COMPARE
457.43967
169.29945
COMPARE
168.28136
62.585
COMPARE
61.90376
23.119616
COMPARE
22.763624
8.4966755
COMPARE
8.348415
OVER
OVER
从这里可以看出结果是可以信赖的,也就是说伯努利方法用在这里也是没有问题的;尽管其中
能用这种解法的条件相当隐晦;
另外在其中我曾经有过一个错误,关于计算
(pow e
(* 2
x))
的值的时候,因为这里x的值有可能比较大,比如20,而采用
(calc 10 (* 2
x))
的方式是计算不出正确结果的;
如果采用以下方式
(defun formula (x)
(sqrt
(+ -0.5
(* c
(calc 40
(* 2
x))))))
又是可以计算出正确结果的;
假设微分方程模型为y'=y+1/y
方程两边分别除以1/y有:
y'/{1/y}=y/{1/y}+1
Go
y'/y^(-1)=1/y^(-2)+1
按照伯努利方法设置v=1/y^(-2),则v'=2*y'*/y^(-1),
GO
v'/2=v+1
Go
v'=2v+2
Go
v'-2v=2
左边部分按照标准解法有:
|u|=C*e^(-2x)
GO
可以得到v=-1/2+C*e^(2x)
这样可以得到|y|=sqrt v =sqrt {-1/2+C*e^(2x)} ,
现在假设初始条件是y(0)=3,那么C=9.5
¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥
下面写程序来证明:
(setq c 9.5)
(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))
y'=y+1/y
(defun expr (x1 y1 )
(+ y1
(* 0.01
(+ y1
(/ 1
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)
(sqrt
(+ -0.5
(* c
(pow e
(* 2
x))))))
(defun exprhelp (x1 )
(if (< (abs (- x1 0))
0.01)
(formula 0)
(expr x1
(exprhelp (- x1
0.01)))))
(defun test (n)
(if (> n 0)
(progn
(print (exprhelp n))
(print 'compare)
(print (formula n))
(test (- n 1)))
(print 'over)))
[25]> (test 20)
1.3757553E9
COMPARE
1.4953814E9
5.08632E8
COMPARE
5.5012006E8
1.8804696E8
COMPARE
2.0237784E8
6.9523024E7
COMPARE
7.445064E7
2.5703448E7
COMPARE
2.7388858E7
9502851.0
COMPARE
1.0075797E7
3513310.5
COMPARE
3706678.5
1298910.3
COMPARE
1363610.8
480221.44
COMPARE
501644.34
177543.2
COMPARE
184544.63
65639.7
COMPARE
67890.164
24267.734
COMPARE
24975.395
8972.051
COMPARE
9187.934
3317.067
COMPARE
3380.0518
1226.3562
COMPARE
1243.4513
453.3967
COMPARE
457.43967
169.29945
COMPARE
168.28136
62.585
COMPARE
61.90376
23.119616
COMPARE
22.763624
8.4966755
COMPARE
8.348415
OVER
OVER
从这里可以看出结果是可以信赖的,也就是说伯努利方法用在这里也是没有问题的;尽管其中
能用这种解法的条件相当隐晦;
另外在其中我曾经有过一个错误,关于计算
(pow e
(* 2
x))
的值的时候,因为这里x的值有可能比较大,比如20,而采用
(calc 10 (* 2
x))
的方式是计算不出正确结果的;
如果采用以下方式
(defun formula (x)
(sqrt
(+ -0.5
(* c
(calc 40
(* 2
x))))))
又是可以计算出正确结果的;