1

我有两个关系：

y=y0+V*t+sin(w*t)    #relation1
dy/dt=(current/P)-(M1/P)*sin(y)+(M2/P)*sin(y+E)+U*cos(y)*sin(w*t)    #relation2

（M1、P、M2、E、w 和 U 是数字常数）我的目标是找到不同电流的 V（电压）。为了做到这一点，我必须为不同的电流数值求解关系 2，并得到 dy/dt，然后使用 y 和 V 之间的关系 <∂y/∂t>=V (<....>表示时间平均值），我必须找到 V。考虑到我不知道 dy/dt 的值。我试过这个

current = 6e-7 : 1e-8 : 8.5e-7;
for k=1:length(current)
f = @(y, t, M1, P, M2, E) (current(k)/P)-(M1/P)*sin(y)+(M2/P)*sin(y+E)+U*cos(y)*sin(w*t);
[t{k}, y{k}] = ode45(f,tspan,y0);
end

这给了我一个细胞中不同电流的y。

我发现下面的代码会给我 dy/dt：

ydot=y(:,2)   #if I use 1 instead of 2 it will give me y)

但是现在，我的问题变成了这样：当我使用这段代码时，它只会给我 1 个电流的 dy/dt，我怎样才能得到不同电流的 dy/dt？

4

2 回答 2

1

试着让你的输出成为一个向量，其中第一个元素是你感兴趣的状态（y），第二个元素是它对时间的导数（dy/dt）；所以y0 = [0;0];或无论你的起始条件是什么。然后为您的 ODE 创建一个单独的文件，我们称之为“myFcn”：

function dydt = myFcn(~, y, M1, P, M2, E, current)      % the ~ is because we are not explicitly dependent on tspan

% Initialize the d/dt vector of our states, y
dydt = zeros(size(y));

% Update the d/dt vector of our states
dydt(1) = y(2);                                         % because (d/dt)y(1) = y(2) = dydt
dydt(2) = current/P - (M1/P)*sin(y) + (M2/P)*sin(y+E);  % your update equation

现在只需替换​​您的句柄并ode45在上面调用：

f = @(y, M1, P, M2, E, current(k))myFunc(y, M1, P, M2, E, current(k));
[t{k}, y{k}] = ode45(f, tspan, y0);

然后，您的输出将是y给出“位置”状态y(1)和“速度”状态的向量y(2)。这些将像以前一样存储在您的单元格数组中。

编辑：

更新了包含current(k)的代码，与 OP 的代码保持一致。

于 2013-09-30T15:24:17.187 回答
-1

如果您想简单区分，可以使用 f2 = jacobian(f,[dt]) 和 matlabfunction(f2) 来取回函数句柄...

于 2013-09-30T15:08:41.713 回答