solving ode in a given grid

I need to solve the system of ode

where t is the parameter. I know the analytical form of the matrix . However the parameter t is the function of x given in the discrete points , i.e. I know only the values . If I use the function in the form ode45(@rhs, [x1, x2, ... xn], initial condition), MATLAB should calculate the rhs of the equation, i.e. the matrix in the points . However, how to let MATLAB know that it should also use the propper values of the parameter t in certain grid points? I mean, I need to have , where .

For example, I solve the equation , where , with . The exact solution is . The main code is

clear all
global t
x = 0:.1:1;
t = x.^2;
[X,Y] = ode45(@test45rhs,x,1);
plot(Y)

I create a function

function dy = test45rhs(x,y)
global t
dy = - t.*y./x;
%dy = - x.*y;

It does not work.

However, if I modify the function

function dy = test45rhs(x,y)
global t
%dy = - t.*y./x;
dy = - x.*y;

Everything works.

[Jan]

What about calling the integration inside a loop over k?

[Koschey Bessmertny]

I guess it is the only solution. Thanks

Answer 1

[Bjorn Gustavsson]

If you have the analytical function for how your t depends on x you should be able to calculate t for any x inside your ODE-function. A modification of your first example along these lines ought to do it:

function dy = test45rhs(x,y)
t = x^2;
dy = - t.*y./x;
%dy = - x.*y;

If you only have t sampled at a set of points x_i then you might do something like this (provided that the underlying relation between t and x is smooth enough):

function dy = test45rhs(x,y,x_i,t_i)
t = interp1(x_i,t_i,x,'pchip');
dy = - t.*y./x;
%dy = - x.*y;

For this case you'll also have to modify the call to ode45 to something like this:

clear all
global t
x = 0:.1:1;
x_i = x;    % for readability
t = x_i.^2; % Just using your illustrating example 
[X,Y] = ode45(@(x,y) test45rhs(x,y,x_i,t_i),x,1);
plot(Y)

Now the "best-interpolated" value of t will be calculated for every call of test45rhs at "the current x-value".

HTH

[Bjorn Gustavsson]

Did this answer solve your problem?

[Koschey Bessmertny]

Thanks for the answer. The interpolation is slow in matlab. I decided to use the Euler integration in a loop.