How to code a couple ODE with nodes

Hello everyone. I am trying to model a thermal system. It is a system of the equation that must be solved:

Xi and Yi are the temperatures at a hypothetical element. (i-1) and (i+1) refer to the previous and next node and A to G are some known coefficients. I consider the ode45 to solve these coupled equations, but the problem is with defining a new parameter such as Z to involve both X and Y, how the nodes (i-1 i i+1) be considered?

Can the following code be run with a few changes?

function zprime = myfun(t,z);

zprime=[A*x(1)(i) + B*x(1)(i-1) + C*x(1)(i+1) + D*x(2)(i) ; E*x(2)(i) + F*x(2)(i-1) + G*x(1)(i)]; % This is absoloutly wrong call!

z0=[a b];

tspan=[0,20];

[t,x]=ode45(@myfun,tspan,z0);

Many thanks!

Answer 1

[Davide Masiello]

The way you describe your problem looks similar to a discretization of a PDE.

You could try to modify the function as follows

function zprime = myfun(t,z);
N = 30;                                             % Number of nodes
A = 1; B = 1; C = 1; D = 1; E = 1; F = 1; G = 1;    % Your known constants

x = z(1:N);                                         % Define variable x as first 30 elements of z
y = z(N+1:2*N);                                     % Define variable y as second 30 elements of z

dxdt(1)       = DEFINE BC;
dxdt(2:end-1) = A*x(2:end-1)+B*x(1:end-2)+C*x(3:end)+D*y(2:end-1);
dxdt(end)     = DEFINE BC;

dydt(1)       = DEFINE BC;
dydt(2:end-1) = E*y(2:end-1)+F*y(1:end-2)+G*x(2:end-1);
dydt(end)     = DEFINE BC;

zprime=[dxdt;dydt];
end

As you can see, the code above is incomplete because you have to define the boundary conditions.

The reason your formulae would not work at, say, i=1 is that i-1=0, which matlab does not accept as index.

Also, at i=30 you have i+1=31 which is larger than the size of x and y.

Hope that clarifies it a bit.

[ehsan rezaei]

Dear Davide, Thank you a lot for your response.

It is a simple code after your correction:

---------------------------------

function zprime = myfun(t,z)

N = 30; % Number of nodes

A = 1; B = 1; C = 1; D = 1; E = 1; F = 1; G = 1; % Your known constants

x = z(1:N); % Define variable x as first 30 elements of z

y = z(N+1:2*N); % Define variable y as second 30 elements of z

dxdt(1) = 2;

dxdt(end) = 5;

dxdt(2:end-1) = A*x(2:end-1)+B*x(1:end-2)+C*x(3:end)+D*y(2:end-1);

dydt(1) = 1;

dydt(2:end-1) = E*y(2:end-1)+F*y(1:end-2)+G*x(2:end-1);

dydt(end) = 3;

zprime=[dxdt;dydt];

end

---------------------------------

z0=[2 1];

tspan=(0:0.1:20);

[t,z]=ode45(@myfun,tspan,z0)

But there is still some error:

Index exceeds the number of array elements (2).

Error in myfun (line 4)

x = z(1:N); % Define variable x as first 30 elements of z

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Untitled6 (line 4)

[t,z]=ode45(@myfun,tspan,z0)

[Davide Masiello]

What is the length of z0?

[ehsan rezaei]

The z0 was selected hypothetically. Its refer to x0 and y0. I want to know where is the problem?

[Davide Masiello]

You need the array of initial condition z0 to be Nx2 long.

You can fix this problem by modifying the definition of N in the function to

N = length(z)/2;

I also noticed something else. You have assigned

dxdt(1)       = 2;
dxdt(end)     = 5;

dydt(1)       = 1;
dydt(end)     = 3;

dxdt and dydt are values of the time derivative of x and y, not the values of x and y themselves. Are you sure they are constants?

If you need to set the value of x and y themselves, than you need to solve algebraic equations at the first and last nodes, and your system switches to a DAE. Matlab can handle this, but you'd need to change ODE solver.

Answer 2

[ehsan rezaei]

Dear Davide,

Thanks for your recommendations, the problem solved!