Длина вектора, возвращаемого функцией COPY_OF_ODEFUN, равна 1, но длина вектора начальных условий равна 2.

Я пытался решить Copy_of_odefun с помощью ode45. И я получил ошибку: длина вектора, возвращаемого 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)))];

Всего 1 ответ


Проблема в глобальных переменных. Вы должны отметить их как в функции, так и в основном скрипте. Просто добавь

global H lambda Cd

на основной сценарий, и он будет работать.

Есть также проблемы со вторым параметром для ode45 - это должно быть просто [8, 100] , а не вектор со всеми необходимыми значениями t . Некоторые хорошие примеры можно найти в документации . В коде могут быть некоторые другие проблемы - он дает NaN для большинства значений y .


Есть идеи?

10000