Matlab: Using interp2 in ode45

Question

I have a 2-dimensional matrix with discrete values and I want to use the ode45 command. To do this I used interp2 to create a function ode45 can use.

The problem is, it looks like ode45 is moving out of my defined area and so interp2 returns NaN-values. To get rid of the NaN I used an extrapolation value but now it seems that ode45 integrates from my initial values to this extrapolation value, disregarding any of my given values in the matrix.

Here is a little example with a 2x2 matrix:

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]

What am I doing wrong? Thanks for any help in advance!

You're passing a constant matrix to ode45 as the first input. I think you mean to call it as [T,C] = ode45(B, tspan, [Xi,Yi]); since B contains your function handle. — Doug Lipinski
When I'm using B instead of @(Xi,Yi)B(Xi,Yi) I'm getting the same results. I'm using Matlab R2010b by the way, don't know if it makes a difference. — user3465431
My mistake, you're actually passing a function that is equivalent to just passing B. I'm actually not sure what 's going on. There's something odd about what your doing though. If B is your ODE and you apparently have two state variables (Xi,Yi) I would think B needs to return 2 values which correspond to the rates of change of Xi and Yi. My best guess is that ode45 is reshaping those arrays in a way that doesn't give an error, but isn't what you want. — Doug Lipinski
Yes, B(Xi,Yi) gives me a new matrix with the interpolated values from matrix A. I don't know what ODE is doing with it, because I thought it would give me the trajectories depending on my values in A. — user3465431

Doug Lipinski Doug Lipinski · Accepted Answer · 2014-03-27T18:03:41

I don't think this will work as you have set it up. Your ODE function B should take inputs of time, and a single input (column) vector of state variables and return a column vector of rates of change for the state variables. Please look at the help documentation for ode45.

If your state variables are coordinate pairs (x,y), you need to return dx/dt and dy/dt for each pair. I'm not sure how this could come from interpolating on a single array. I think you need something like the following:

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

Obviously this code relies heavily on managing array shapes to maintain the correct ordering for x, y, dx/dt, and dy/dt. There may be a simpler way to write it.

EDIT
The key here is to clearly define a state vector and define your ODE function to match that state vector. The state vector must be a column vector and your ODE function must return a column vector. In the code above I have chosen to represent the state of the system as a vector formated as states(:,1) = [x1 y1 x2 y2 x3 y3 ... ]. That means your ODE must return a column vector of the form

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

You will also require 2 interpolations, 1 for the x components and 1 for the y, to get the rates of change based on Ax and Ay. The way I chose to do this was with the line

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));

This line is a bit complex and difficult to understand because it is written as an anonymous function. If you define a stand-alone function for this it will be much more clear and could be written as

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);

One last note:
To use ode45, your ODE function must take inputs of a time (t above) and a state vector, even if you don't use the time. Additional arguments are optional, see the documentation on ode45 for more details.

Matlab: Using interp2 in ode45

1 Answers