Matlab：在ode45中使用interp2

Question

我有一個具有離散值的2維矩陣，我想使用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必須交回

​​

您還需要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上的文檔。