Method ode45 applied to set of n 1st order ODE

Hi everybody. Not new to MatLab but facing a problem.

Let's say I have the following set of n=5 ODEs (not interested in initial/boundary conditions, as I know them):

1) dx1/dt = f1(x1(t),x2(t)) = c1*(x2-x1);
2) dx2/dt = f2(x1(t),x2(t),x3(t)) = c2*(x1+x3-2*x2);
3) dx3/dt = f3(x2(t),x3(t),x4(t)) = c3*(x2+x4-2*x3);
4) dx4/dt = f4(x3(t),x4(t),x5(t)) = c4*(x3+x5-2*x4);
5) dx5/dt = f5(x4(t),x5(t)) = c5*(x4-x5);

Let's assume the initial conditions are expressed by the vector (array):

x0 = [x1(0),x2(0),x3(0),x4(0),x5(0)].

I can solve this with the following code:

f=@(t,x) [c1*(x2-x1);c2*(x1+x3-2*x2);c3*(x2+x4-2*x3);c4*(x3+x5-2*x4);c5*(x4-x5)] %%%c1,..,c5 are constants
[t,xsol] = ode45(f,[0,tend],x0); %%%[0,tend] is the range in which I want to solve;

So far so good. Problem arises when I don't write explicitily the equations in brackets in the line that defines the anonymous function f.

Let's say the user is required to input at runtime the number of equations n, and that they all have the same aspect of the ones written above. Of course I could rewrite the equations from scratch, but what if n is very big, e.g. n =100? And if n is frequentl changed?

Using a for loop I was able to define the vector of coefficient CC = [c1,c2,...cn] for any given integer n.

Inside the same for loop I was able to define the vector XX = [x2-x1;x1+x3-2*x2;x2+x4-2*x3;...;x(n-1)-xn]

(Sorry that I had to put parentheses for the last two elements, but xn-1 would be less understable, I guess).

Of course, above the loop I declared x1,...xn as symbolic variable as follows:

sym('x'[1,n])

It seems to work well, except that to accumulating values from the for loop I could preallocate the size of the coefficient vector CC (declaring a zeros(n,1)) but I wasn't able to do the same with the XX vector (and I had to define it as XX = [] or, at most, as cell array XX = cell(n,1)).

Anyway, then I compute the element-wise multiplication between CC and XX, in order to get the column vector that is needed in the anonymous function, and I get it correct but it's of sym class.

argument = XX.*CC; %%%sorry if I am making mistakes with transposing here but it's not the point;
f = @(t,x) argument;
[t,xsol] = ode45(f,[0,tend],x0); %%%here x0 = [x1,...xn] of course

This doesn't work anymore and MatLab complains about odefunction arguments being not double or single. I am at a dead end, here. Perhaps I should try to substitute the anonymous function with a function defined in a separate m-file, but I have no clue of how to do it.

Could you please help me? It's very important. As a side question, I would like to know why if I get argument as a row vector (argument is symbolic) and then I transpose it the usual way (' mark), I get conj(x1), conj(x2), instead of the same vector in column form.

Thank you very much to anybody who will spend time at reading and trying to help. Sorry for mistakes related to English language, but I am not a native speaker and never had the chance to spend some time abroad.

Nicola

Answer 1

[Mischa Kim]

Hi Nicola, use the below to get started. Essentially you want to define the differential equations symbolically and then convert the system to solve numerically.

• The matrix M below is not quite correct. Please take some time to double-check.
• The final command creates the function handle you can use to solve the differential equation, e.g. with ode45.

function fh = flex_ode(numii)
syms t

% create symbolic variables for x and c
x   = sym('x',[1 numii]);
c   = sym('c',[1 numii]);

% create RHS of differential equation
M   = diag(ones(numii-1,1),1) + diag(ones(numii-1,1),-1) - 2*diag(ones(numii,1));
RHS = (c.').*(M*x.');

% create symbolic functions, x(t), and their derivatives
xt  = arrayfun(@(ii) symfun(['x' int2str(ii) '(t)'],t),1:numii,'UniformOutput',false);
dxt = arrayfun(@(ii) diff(xt{1,ii},t),1:numii,'UniformOutput',false);

RHS = subs(RHS,x,xt);

% define differential equation
DE  = [];
for ii = 1:numii
    DE = [DE; dxt{ii} == RHS(ii)]; %#ok<AGROW>
end

% compute vector field and MATLAB function
[VF,Ysubs] = odeToVectorField(DE);
fh         = matlabFunction(VF,'vars',[{'t','Y'}, {c.'}],'outputs', {'dYdt'});
end