Matlab:在 ode45 中使用 interp2

Question

我有一个带有离散值的二维矩阵，我想使用 ode45 命令。 为此，我使用 interp2 创建了一个 ode45 可以使用的函数。

问题是，看起来 ode45 正在移出我定义的区域，因此 interp2 返回 NaN 值。为了摆脱 NaN，我使用了外推值，但现在看来 ode45 从我的初始值积分到这个外推值，忽略了矩阵中我给定的任何值。

这是一个 2x2 矩阵的小例子:

clear all;
close all;
clc;

A = rand(2, 2);              % the matrix with my values
x = 1:2;         
y = x;
[X, Y] = meshgrid(x, y);     % grid on which my values are defined
xi = 1:0.5:2;              
yi = xi;
[Xi, Yi] = meshgrid(xi, yi); % interpolation grid

tspan = 0:0.2:1;

B = @(xq,yq)interp2(X, Y, A, xq, yq, 'linear', 10);   % interpolation function with extrapval of 10

[T,C] = ode45(@(Xi,Yi)B(Xi,Yi), tspan, [Xi,Yi]);       % ode45 function using interp-function and intital values [Xi,Yi]

我做错了什么？ 感谢您提前提供的任何帮助!

Answer 1

我认为这不会像您设置的那样起作用。您的 ODE 函数 B 应接受时间输入和状态变量的单个输入(列)向量，并返回状态变量变化率的列向量。请查看ode45的帮助文档。

如果您的状态变量是坐标对 (x,y)，则需要为每对返回 dx/dt 和 dy/dt。我不确定这是如何通过对单个数组进行插值得出的。我认为您需要类似以下内容:

clear

%// Define the right hand side of your ODE such that
%// dx/dt = Ax(x,y)
%// dy/dt = Ay(x,y)
%// For demonstration I'm using a circular velocity field.
xref = linspace(-pi,pi);
yref = linspace(-pi,pi);
[xref yref] = meshgrid(xref,yref);

r=sqrt(xref.^2+yref.^2);
Ax = [-yref./r];
Ay = [xref./r];

%// Initial states:
xi = 1:0.05:2;
yi = xi;
[Xi, Yi] = meshgrid(xi, yi); % interpolation grid

%// state array = [x1; y1; x2; y2; ... ]
states = [Xi(:)'; Yi(:)'];
states = states(:);

B = @(t,states) reshape([ interp2(xref,yref,Ax,states(1:2:end),states(2:2:end)) ...
                          interp2(xref,yref,Ay,states(1:2:end),states(2:2:end))]',size(states));

tspan=0:.02:10;
[T,C] = ode45(B, tspan, states);

for i=1:size(C,1)
    figure(1)
    clf
    plot(C(i,1:2:end),C(i,2:2:end),'k.')
    axis(pi*[-1 1 -1 1])
    drawnow
end

显然，这段代码在很大程度上依赖于管理数组形状来维护 x、y、dx/dt 和 dy/dt 的正确顺序。可能还有更简单的写法。

编辑
这里的关键是明确定义状态向量并定义 ODE 函数以匹配该状态向量。状态向量必须是列向量，并且 ODE 函数必须返回列向量。在上面的代码中，我选择将系统状态表示为格式为 states(:,1) = [x1 y1 x2 y2 x3 y3 ... ] 的向量。这意味着您的 ODE 必须返回以下形式的列向量

[ d/dt(x1) ]
[ d/dt(y1) ]
[ d/dt(x2) ]
[ d/dt(y2) ]
[ d/dt(x2) ]
[ d/dt(y2) ]
[   ...    ]
[   ...    ]

您还需要进行 2 次插值，1 次用于 x 分量，1 次用于 y 分量，以获得基于 Ax 和 Ay 的变化率。我选择这样做的方式是使用该行

B = @(t,states) reshape([ interp2(xref,yref,Ax,states(1:2:end),states(2:2:end)) ...
                          interp2(xref,yref,Ay,states(1:2:end),states(2:2:end))]',size(states));

这行代码有点复杂且难以理解，因为它是作为匿名函数编写的。如果你为此定义一个独立的函数，它会更清晰，可以写成

function ODEvals = B(t,states,xref,yref,Ax,Ay)
x(:,1) = states(1:2:end); %// extract x values from states as a column vector
y(:,1) = states(2:2:end); %// extract y values
dxdt(:,1) = interp2(xref,yref,Ax,x,y); %// interpolate Ax
dydt(:,1) = interp2(xref,yref,Ay,x,y); %// interpolate Ay

%// concatenate the results, dxdt and dydt are column vectors
%//  1) put the as column 1 and 2
%//  2) take the transpose so they become rows one and two:
%//         [d/dt(x1)  d/dt(x2) ... ]
%//         [d/dt(y1)  d/dt(y2) ... ]
%//  3) reshape into a single column, the ordering will be:
%//        [ d/dt(x1) ]
%//        [ d/dt(y1) ]
%//        [ d/dt(x2) ]
%//        [ d/dt(y2) ]
%//        [ d/dt(x2) ]
%//        [ d/dt(y2) ]
%//        [   ...    ]
%//        [   ...    ]
ODEvals = reshape([dxdt dydt]',[],1);

最后一点:
要使用 ode45，您的 ODE 函数必须接受时间(上面的 t)和状态向量的输入，即使您不使用时间也是如此。其他参数是可选的，请参阅 ode45 上的文档了解更多详细信息。