COPY_OF_ODEFUN返回的向量长度为​​1，但初始条件向量的长度为2

Question

我尝试通过 ode45 解决 Copy_of_odefun。我得到了错误：COPY_OF_odefun返回的向量长度为​​1，但初始条件向量的长度为2。我的代码如下。

主要功能：

clear all; clc;
%% Inpput parameters
lambda = 0.0625;
H = 2.3;
zw_H = 2;
ztop_H = 8;
Cdh = 1.2;
a = 9.6 * lambda;
Cd = Cdh * (1 - exp(-2 * a));
d_H = 0.066;  %% Macdonald, 2000(lambda = 0.05)
%% Running
t = linspace(8,0.01,100);
y0 = [1.75; 0.05];
[t,y] = ode45(@Copy_of_odefun,t,y0);

copy_of_odefun 函数：

function dy = Copy_of_odefun(t,y)
global H lambda Cd
%% mixing length closure parameter
% within canopy
z_H1 = [0.01 1];
lc_H = [0.2463 0.2463];
lc_dz = [0 0];
% shear layer
z_H2 = [1.01 2];
lm_H = [0.2516 0.7736];
lm_dz = [0.5273 0.5273];
% inertial layer
z_H3 = [2.01 8];
lz_H = [0.7776 3.1736];
lz_dz = [0.4 0.4];
% for whole computation domain
zm_H = [z_H1 z_H2 z_H3];
l3_H = [lc_H lm_H lz_H];
l3_dz = [lc_dz lm_dz lz_dz];
% interpolate to obtain value at each point of t
lx_H = interp1(zm_H,l3_H,t);
lx_dz = interp1(zm_H,l3_dz,t);

% second order translated into one order
dy = [y(2); (-lx_H * H * lx_dz * y(2)^2 / (H^2 * lx_H^2 * y(2)) + 0.25 * y(1)^2 * lambda / H * Cd / (H^2 * lx_H^2 * y(2)))];

Answer 1

问题出在全局变量中。您需要在函数和主脚本中标记它们。只需添加

global H lambda Cd

到主脚本，它将运行。

ode45 的第二个参数也有问题 - 它应该只是 [8, 100]，而不是所有需要的值 t 的向量。在documentation 中可以找到一些很好的例子。代码中可能还有其他一些问题 - 它为大多数 y 值提供了 NaN。