interp2在ode45中的应用

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) ]
[   ...    ]
[   ...    ]

您还需要两个插值，一个用于x组件，一个用于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的文档。