【微分方程】最佳砍伐强度问题

KeyTo9_Fans

一开始有10000吨硬木树和6000吨软木树，每过dt年，都会新长出(0.1*dt*H)吨硬木树和(0.25*dt*S)吨软木树
而在这期间，会有(0.00001*dt*H^2)吨硬木树和(0.00004*dt*S^2)吨软木树因种群内竞争而灭亡
此外，还另有(0.000005*dt*H*S)吨硬木树和(0.00002*dt*S*H)吨软木树因种群间竞争而灭亡
其中，dt是一个趋于0的无穷小量，H是硬木树当前时刻的总数量（吨），S是软木树当前时刻的总数量（吨）

问题1：假设每过dt年，都会有dt*x吨硬木树和dt*y吨软木树被砍伐掉，其中x和y是砍伐强度，单位是（吨/年）。
问当x和y满足什么条件时（满足该条件的区域称为可行域），永远都不会有任何一种树被砍光？请在x-y平面坐标系里画出该可行域。

问题2：假设硬木树的单价是软木树的k倍，问当k=2时，x和y应分别取何值，才能使得永远都不会有任何一种树被砍光且年利润最大？

问题3：当k发生变化时，问题2的答案x(k)和y(k)是如何随k变化的？

Jack315

微分方程：
\(\frac{d}{dt}h(t)=a_1h(t)-a_2h^2(t)-a_3h(t)s(t)-x\)
\(\frac{d}{dt}s(t)=b_1s(t)-b_2s^2(t)-b_3h(t)s(t)-y\)
其中：
\(h(t)\) 是时刻 \(t\) 时的硬木树总量；\(s(t)\) 是时刻 \(t\) 时的软木树总量。
\(x\) 为硬木树的砍伐强度；\(y\) 为软木树的砍伐强度。
\(\begin{matrix}
a_1=0.1&a_2=0.00001&a_3=0.000005\\
b_1=0.25&b_2=0.00004&b_3=0.00002
\end{matrix}\)
初始条件：
\(\begin{matrix}h(t)\vert_{t=0}=10000&s(t)\vert_{t=0}=6000\end{matrix}\)

Simulink 模型文件：
砍伐强度_Simulink.rar (22.51 KB, 下载次数: 0)

2024-12-27 17:57 上传

点击文件名下载附件

硬木树被砍光，软木树疯长 \((x=165,y=20)\) ：

软木树被砍光，硬木树疯长 \((x=90,y=45)\) ：

KeyTo9_Fans

怎么会疯长呢？不是有(0.00001*dt*H^2)吨硬木树和(0.00004*dt*S^2)吨软木树因种群内竞争而灭亡吗？达到环境容纳量它就增长不了了

Jack315

KeyTo9_Fans 发表于 2024-12-27 20:26
怎么会疯长呢？不是有(0.00001*dt*H^2)吨硬木树和(0.00004*dt*S^2)吨软木树因种群内竞争而灭亡吗？达到环境 ...

是模型设置问题：
在模型内，一种树砍光之后不是维持在零，而是成了负数。设置下限值为 0 应该就行了。
【更正图】
\((x=165,y=20)\)硬木树被砍光，软木树疯长至环境容纳量：

\((x=90,y=45)\)软木树被砍光，硬木树疯长至环境容纳量：


可行域：

兰色区域是硬木树被砍光；红色区域是软木树被砍光。
在可行域的尖角处（~180，~90）右上角，【可能】是一条脆弱的平衡曲线（不稳定不动点）。

计算可行域的 Matlab 代码：

1. clear
   
2. w = warning;
   
3. warning('off');
   
4. 
   
5. nh = 200;
   
6. ns = 100;
   
7. hdata = zeros(nh, ns, 3);
   
8. sdata = zeros(nh, ns, 3);
   
9. tspan = [0 1000];
   
10. y0 = [10000 6000];
    
11. 
    
12. for h = 1 : nh
    
13.     for s = 1 : ns
    
14.         [t,y] = ode45(@(t,y) odefcn(t,y,h-1,s-1), tspan, y0);
    
15.         hdata(h, s, 1) = h;
    
16.         hdata(h, s, 2) = s;
    
17.         hdata(h, s, 3) = CrossZero(t, y(:,1));
    
18. 
    
19.         sdata(h, s, 1) = h;
    
20.         sdata(h, s, 2) = s;
    
21.         sdata(h, s, 3) = CrossZero(t, y(:,2));
    
22.     end
    
23. end
    
24. 
    
25. warning(w);

复制代码

微分方程函数：

1. function dydt = odefcn(t,y,h,s)
   
2.     dydt = zeros(2,1);
   
3.     dydt(1) = 0.1*y(1)-0.00001*y(1)^2-0.000005*y(1)*y(2)-h;
   
4.     dydt(2) = 0.25*y(2)-0.00004*y(2)^2-0.00002*y(1)*y(2)-s;
   
5. end

复制代码

过零点数据提取函数：

1. function tzero = CrossZero(t, y)
   
2.     ymin = min(y);
   
3.     if ymin < 0
   
4.         k = find(y > 0, 1, 'last');
   
5.         % 过零时间。
   
6.         t0 = t(k); t1 = t(k + 1);
   
7.         y0 = y(k); y1 = y(k + 1);
   
8.         tzero = (t1 - t0) / (y1 - y0) * (0 - y0) + t0;
   
9.     else
   
10.         tzero = nan;
    
11.     end
    
12. end

复制代码

【猜想】
在特定的模型参数下，有可能导致 \(h(t)\) 和 \(s(t)\) 出现振荡，形成二维“吸引子”。

Jack315

可行区域由两条边界曲线和两坐标轴所围成，如下图所示：

边界曲线方程为多项式：
\(y=\sum_{i=0}^n\beta_ix^i\)，硬木树 \(n=2\)，软木树 \(n=5\) 。
其系数如下表所示：
\(\begin{array}{|c|c|c|}
\hline \text{项}&\text{硬木树}&\text{软木树}\\
\hline \text{截距}&-4040.6434468525&22.7007398066704\\
\hdashline x&43.5778508771939&0.0773407651352562\\
\hdashline x^2&-0.11448658410733&0.00447486018733332\\
\hdashline x^3&-&-0.0000601214112270588\\
\hdashline x^4&-&4.18420901158062E-07\\
\hdashline x^5&-&-9.46834887323358E-10\\
\hline \end{array}\)
边界曲线四个交点坐标如下表所示：
\(\begin{array}{|c|c|c|}
\hline \text 交点&x&y\\
\hline O&0&0\\
\hdashline A&159.911&0\\
\hdashline B&0&22.7007\\
\hdashline C&178.78&90.533\\
\hline \end{array}\)

Jack315

问题 2 和 3 审题被卡住了。
这是一个有约束条件的优化问题，但目标函数不知道如何建立？
有没有高手可以给出这个目标函数？

KeyTo9_Fans

我按dt=0.001（年）模拟出来的尖角坐标是(174.212, 86.888)，dt继续缩小要耗费我大量的cpu（如有新进展，我会更新此楼层）

如果不考虑非稳定平衡状态，那么第2、3问的最优解很可能就是在这个尖角上（无论k取多少）

Jack315

5# 的边界曲线是用有误差的数据拟合出来的，所以是缺乏精确度的。

只要价格与成交量（砍伐强度 x, y）无关，最优解都在这个尖角处，而与 k 无关。

KeyTo9_Fans

我用MATLAB求得的尖角上的稳定状态如下：

H = 5107.5072182050216775844458966333
S = 2963.178504298315242025548166767
x = 174.21214398267014760652871941207
y = 86.888440150482685819608142236456

代码和运行结果如下：

1. %MATLAB代码：
   
2. >> syms H S x y
   
3. >> [H S]=solve(0.1*H-0.00001*H^2-0.000005*H*S-x==0,0.25*S-0.00004*S^2-0.00002*S*H-y==0,H,S)
   
4. >> vpa([subs(H,[x,y],[174.213,86.889]);subs(S,[x,y],[174.213,86.889])])
   
5. >> H=real(ans(3)),S=real(ans(7))
   
6. >> x=0.1*H-0.00001*H^2-0.000005*H*S,y=0.25*S-0.00004*S^2-0.00002*S*H
   
7. %上述代码的运行结果如下：
   
8. H =
   
9. (55340232221128653339*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 1)^3 - 438110171750601841448750*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 1)^2 - 2305843009213693952000000000*y + 1844674407370955161600000*x*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 1) + 1383505805528216333475000*y*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 1) + 576460752303423488000000000*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 1))/(922337203685477580800000*y)
   
10. (55340232221128653339*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 2)^3 - 438110171750601841448750*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 2)^2 - 2305843009213693952000000000*y + 1844674407370955161600000*x*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 2) + 1383505805528216333475000*y*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 2) + 576460752303423488000000000*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 2))/(922337203685477580800000*y)
    
11. (55340232221128653339*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 3)^3 - 438110171750601841448750*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 3)^2 - 2305843009213693952000000000*y + 1844674407370955161600000*x*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 3) + 1383505805528216333475000*y*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 3) + 576460752303423488000000000*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 3))/(922337203685477580800000*y)
    
12. (55340232221128653339*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 4)^3 - 438110171750601841448750*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 4)^2 - 2305843009213693952000000000*y + 1844674407370955161600000*x*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 4) + 1383505805528216333475000*y*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 4) + 576460752303423488000000000*root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 4))/(922337203685477580800000*y)
    
13. S =
    
14. root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 1)
    
15. root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 2)
    
16. root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 3)
    
17. root(z^4 - (438110171750601841448750*z^3)/55340232221128653339 + z^2*((1844674407370955161600000*x)/55340232221128653339 + (3228180212899171495075000*y)/55340232221128653339 + 576460752303423488000000000/55340232221128653339) - (4611686018427387904000000000*y*z)/18446744073709551113 + (46116860184273879040000000000*y^2)/55340232221128653339, z, 4)
    
18. ans =
    
19.                                       2349.6575722908332712148514023534
    
20.                                       6601.9946579657899493906413443453
    
21. 5107.5072182050216775844458966333 + 11.321684669266507062525010099168i
    
22. 5107.5072182050216775844458966333 - 11.321684669266507062525010099168i
    
23.                                       471.88579512390400017454358199025
    
24.                                       1518.4238629461322278871678144932
    
25.    2963.178504298315242025548166767 - 7.521638697686384185979495370697i
    
26.    2963.178504298315242025548166767 + 7.521638697686384185979495370697i
    
27. H =
    
28. 5107.5072182050216775844458966333
    
29. S =
    
30. 2963.178504298315242025548166767
    
31. x =
    
32. 174.21214398267014760652871941207
    
33. y =
    
34. 86.888440150482685819608142236456

复制代码

上面这段代码的注释如下：

>> syms H S x y
%这一行代码是把变量 H、S、x、y 定义成符号

>> [H S]=solve(0.1*H-0.00001*H^2-0.000005*H*S-x==0,0.25*S-0.00004*S^2-0.00002*S*H-y==0,H,S)
%假设x和y已知，求解H和S，使得两种树的增长率均为0，解出来的H和S均为一元四次方程的根，方程的系数里均含有x和y

>> vpa([subs(H,[x,y],[174.213,86.889]);subs(S,[x,y],[174.213,86.889])])
%把x和y赋值为尖角外面的值，代入一元四次方程的系数，然后求解H和S，求解结果均为2个实数和2个虚数

>> H=real(ans(3)),S=real(ans(7))
%把H和S赋值为虚数解的实数部分，得到尖角处两种树稳定时的数量

>> x=0.1*H-0.00001*H^2-0.000005*H*S,y=0.25*S-0.00004*S^2-0.00002*S*H
%把尖角处两种树稳定时的数量代入增长率方程，得到尖角处的砍伐强度x和y

上面这段代码是基于下面的观察结果写出来的：

1. 我发现如果x和y取的是尖角内的值(174.212,86.888)，那么解出来的稳定点(H,S)均为实数：
>> vpa([subs(H,[x,y],[174.212,86.888]);subs(S,[x,y],[174.212,86.888])])
ans =
2349.6361404946074689153286260937
6602.0288286952797969994565917793
  5109.320703458368442169474332709
5105.6809940184108676901249893832
471.87851997457151404263183568529
1518.4131308980611019381481630018
2961.9784793205720680257482166775
2964.3965364734620281062795146528

上面这个结果说明如果x和y落在可行域内，我们能找到实数的稳定点(H,S)

2. 如果x和y取的是尖角外的值(174.213,86.889)，那么(H,S)后面的两组解就跑到虚数上去了：
>> vpa([subs(H,[x,y],[174.213,86.889]);subs(S,[x,y],[174.213,86.889])])
ans =
                                      2349.6575722908332712148514023534
                                      6601.9946579657899493906413443453
5107.5072182050216775844458966333 + 11.321684669266507062525010099168i
5107.5072182050216775844458966333 - 11.321684669266507062525010099168i
                                      471.88579512390400017454358199025
                                      1518.4238629461322278871678144932
   2963.178504298315242025548166767 - 7.521638697686384185979495370697i
   2963.178504298315242025548166767 + 7.521638697686384185979495370697i

我的理解是：

当x和y移动到可行域的边界上时，两个靠得很近的实数解就会合并成重根

当x和y移动到可行域外面时，这个稳定点(H,S)就飞到虚数上去了

也就是当实数范围内找不到相应的稳定状态(H,S)来对应砍伐强度(x,y)时，这个(x,y)就不在可行域里

因此我猜测重根H=5107、S=2963的值就是夹角处的两种树稳定后的数量

因此把这个重根稳态(H,S)代入增长率方程后，求得的x和y值就是尖角处的砍伐强度值

不知道我理解得对不对呢？

Jack315

关注最后的稳态不失为一个好的视角，比起动态过程的视角更适合这个题。
方程解的行为值得仔细考察，尤其是在可行域的边界处。
另外方程有四个解，这些解是什么含意也值得探究。