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

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)

6 天前 上传

点击文件名下载附件

硬木树被砍光，软木树疯长 \((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}\)