ode45 with a array of vector

I am solving 6 ODEs simultaneously. My eqns are:

eq1 = diff(x,t) == u1+vp.*p1 ;

eq2 = diff(y,t) == u2+vp.*p2 ;

eq3 = diff(z,t) == u3+vp.*p3 ;

eq1 = diff(p1,t) == a1*p1+a2*p2+a3*p3 ;

eq2 = diff(p2,t) == b1*p1+b2*p2+b3*p3 ;

eq3 = diff(p3,t) == c1*p1+c2*p2+c3*p3 ;

Here, a1, a2, a3, b1, b2 b3, c1, c2, c3 and vp are constants.

u1, u2, and u3 are vector of dimension [1 X 100]. Each values of u1, u2, and u3 corresponds to time points in tValues = linspace(0,10,100). I want to compute p1, p2, p3 and x, y, z for tValues = linspace(0,10,100).

My code is following: ;

vars = [x(t); z(t); y(t); p1(t); p2(t); p3(t)];

V = odeToVectorField([eq2 eq1 eq3 eq4 eq5 eq6]);

M = matlabFunction(V,'vars', {'t','Y'});

y0=[0 0 0 p1_0 p2_0 p3_0];

ySol_a = ode45(M,interval,y0);

Once the code is run, it shows following error message:

MuPAD error: Error: Cannot convert the initial value problem to an equivalent dynamical system. Either the differential equations cannot be solved for the highest derivatives or inappropriate initial conditions were specified. [numeric::ode2vectorfield]

Answer 1

[Torsten]

a1 = ...;
a2 = ...;
a3 = ...;
b1 = ...;
b2 = ...;
b3 = ...;
c1 = ...;
c2 = ...;
c3 = ...;
vp = ...;
tValues = linspace(0,10,100);
u1Values = ...;
u2Values = ...;
u3Values = ...;
u1fun = @(t)interp1(tValues,u1Values,t);
u2fun = @(t)interp1(tValues,u2Values,t);
u3fun = @(t)interp1(tValues,u3Values,t);
fun = @(t,y)[u1fun(t)+vp*y(4);u2fun(t)+vp*y(5);u3fun(t)+vp*y(6);a1*y(4)+a2*y(5)+a3*y(6);b1*y(4)+b2*y(5)+b3*y(6);c1*y(4)+c2*y(5)+c3*y(6)];
tspan = tValues;
p1_0 = ...;
p2_0 = ...;
p3_0 = ...;
y0 = [0 0 0 p1_0 p2_0 p3_0]; 
[T,Y] = ode45(fun,tspan,y0);
plot(T,Y)

[Torsten]

Am I making any mistake?

I don't know. But why don't you solve the differential equations all together instead of calling ode45 two times ?

[Torsten]

Could you please suggest if the above code is technically correct (for Line number 229 onwards) ?

As the code executes without errors, it's technically correct. I cannot say whether your equations are correct.

[Torsten]

Since you have chosen the output times in your second call to ode45 equal to the evaluation times using "deval" of your first call to ode45, u1Values and p1_active correspond to the same time vector, namely "tValues". But I repeat: you shouldn't use two calls to ode45 with interpolation of the results of the first call in the second call. You should integrate the complete set of equations in one call.

[Torsten]

Is there any specific reason for this suggestion?

Sure. Interpolation means loss of precision.

[Torsten]

How do I provide this dynamic initial condition?

I think I answered this already. An initial condition is not dynamic.

If you want to set

v1(t) = constant1 + vp*p1, 
v2(t) = constant2 + vp*p2, 
v3(t) = constant3 + vp*p3

your equations to integrate become

diff(x,t) = constant1 + vp*p1
diff(y,t) = constant2 + vp*p2
diff(z,t) = constant3 + vp*p3
diff(p1,t) == JTx_1+JTx_2+JTx_3;
diff(p2,t) == JTy_1+JTy_2+JTy_3;
diff(p3,t) == JTz_1+JTz_2+JTz_3;

[AkB]

Thanks, it works.