Error on solving differential equation using dsolve

Hi, I have 4 differential equation and I get an error when I use dsolve to solve it.

syms x1(t) x2(t) x3(t) x4(t)
x0=[ 0; 0; 0; 0];
q1=x1(t);
dq1=x2(t);
q2=x3(t);
dq2=x4(t);
dq = [dq1;dq2];
alfa=6.7;beta=3.4;
epc=3.0;eta=0;
m1=1;l1=1;
lc1=1/2;I1=1/12;
g=9.8;
e1=m1*l1*lc1-I1-m1*l1^2;
e2=g/l1;
% global u;
H=[alfa+2*epc*cos(q2)+2*eta*sin(q2),beta+epc*cos(q2)+eta*sin(q2);
 beta+epc*cos(q2)+eta*sin(q2),beta];
C=[(-2*epc*sin(q2)+2*eta*cos(q2))*dq2,(-epc*sin(q2)+eta*cos(q2))*dq2;
 (epc*sin(q2)-eta*cos(q2))*dq1,0];
G=[epc*e2*cos(q1+q2)+eta*e2*sin(q1+q2)+(alfa-beta+e1)*e2*cos(q1);
 epc*e2*cos(q1+q2)+eta*e2*sin(q1+q2)];

 tol(1)=u(1);
 tol(2)=u(2);

 ddq=inv(H)*(tol'-C*dq-G);

 ode1 = diff(x1) == x2(t);
 ode2 = diff(x2) == ddq(1);
 ode3 = diff(x3) == x4(t);
ode4 = diff(x4) == ddq(2);
odes=[ode1;ode2;ode3;ode4];

Now If I just dsolve(odes), it return empty sym. If I use dsolve(odes,x0) it give me an error.

 >> dsolve(odes)
 Warning: Explicit solution could not be found. 
 > In dsolve (line 201) 

ans =

[ empty sym ]

OR

 >> dsolve(odes,x0)
 Error using mupadengine/feval (line 163)
 The equations are invalid.

 Error in dsolve>mupadDsolve (line 336)
 T = feval(symengine,'symobj::dsolve',sys,x,options);

 Error in dsolve (line 193)
 sol = mupadDsolve(args, options);

Answer 1

[Star Strider]

There is probably not an analytic solution to your system.

The best way to proceed is to use matlabFunction to create an anonymous function for your system, then use ode45 or one of the others to integrate it.

Your slightly revised code becomes:

syms x1 x2 x3 x4 t u1 u2 real
% syms x1(t) x2(t) x3(t) x4(t) u1 u2 real
x0=[ 0; 0; 0; 0];
% q1=x1(t);
% dq1=x2(t);
% q2=x3(t);
% dq2=x4(t);
q1 = x1;
dq1 = x2;
q2 = x3;
dq2 = x4;
dq = [dq1;dq2];
alfa=6.7;beta=3.4;
epc=3.0;eta=0;
m1=1;l1=1;
lc1=1/2;I1=1/12;
g=9.8;
e1=m1*l1*lc1-I1-m1*l1^2;
e2=g/l1;

H=[alfa+2*epc*cos(q2)+2*eta*sin(q2),beta+epc*cos(q2)+eta*sin(q2);
 beta+epc*cos(q2)+eta*sin(q2),beta];
C=[(-2*epc*sin(q2)+2*eta*cos(q2))*dq2,(-epc*sin(q2)+eta*cos(q2))*dq2;
 (epc*sin(q2)-eta*cos(q2))*dq1,0];
G=[epc*e2*cos(q1+q2)+eta*e2*sin(q1+q2)+(alfa-beta+e1)*e2*cos(q1);
 epc*e2*cos(q1+q2)+eta*e2*sin(q1+q2)];

   tol(1)=u1;
   tol(2)=u2;

   ddq=inv(H)*(tol'-C*dq-G);

%  ode1 = diff(x1) == x2(t);
%  ode2 = diff(x2) == ddq(1);
%  ode3 = diff(x3) == x4(t);
% ode4 = diff(x4) == ddq(2);
ode1 = x2;
ode2 = ddq(1);
ode3 = x4;
ode4 = ddq(2);
odes=[ode1;ode2;ode3;ode4];

ODESYS = matlabFunction(odes, 'Vars',{t [x1 x2 x3 x4 u1 u2]});

Then use ‘ODESYS’ (call it whatever you wish) with the numeric solver. Here, the column vector ‘in2’ contains all the arguments in 'Vars', corresponding to that variable list. See the documentation for matlabFunction to understand what it does.

[Mary]

Thanks. What is the x0 when I use ode45 to solve it. I tried x0 of length 6 but it does not work.

[Star Strider]

My pleasure.

The initial conditions vector would be length 4.

I’m not sure what ‘u1’ and ‘u2’ are (my guess is that they’re inputs of some description), so I leave it to you to decide how you want to deal with them in the code, and with respect to your function.

[Mary]

Ok, thanks. But it even does not work with x0 of length 4.

     ode45(ODESYS,[0 5],x0);
Index exceeds matrix dimensions.

 Error in
 symengine>@(t,in2)[in2(:,2);(in2(:,5).*
 (-1.7e2./3.0)+cos(in2(:,1)+in2(:,3)).*1.666e3+cos(in2(:,1)).*1.5
 08655555555556e3-in2(:,4).^2.*sin(in2(:,3)).*1.7e2-
 in2(:,2).*in2(:,4).*sin(in2(:,3)).*3.4e2)./(cos(in2(:,3)).^2.*1.
 5e2-1.87e2)-((cos(in2(:,3)).*1.5e1+1.7e1).*(-
 in2(:,6)+cos(in2(:,1)+in2(:,3)).*
 (1.47e2./5.0)+in2(:,2).^2.*sin(in2(:,3)).*3.0).*
 (1.0e1./3.0))./(cos(in2(:,3)).^2.*1.5e2-1.87e2);in2(:,4);
 ((cos(in2(:,3)).*6.0e1+6.7e1).*(-in2(:,6)+cos(in2(:,1)+in2(:,3)).*
 (1.47e2./5.0)+in2(:,2).^2.*sin(in2(:,3)).*3.0).*
 (5.0./3.0))./(cos(in2(:,3)).^2.*1.5e2-1.87e2)+
 ((cos(in2(:,3)).*1.5e1+1.7e1).*(in2(:,5)-
 cos(in2(:,1)+in2(:,3)).*(1.47e2./5.0)-
 cos(in2(:,1)).*2.662333333333333e1+in2(:,4).^2.*sin(in2(:,3)).*3
 .0+in2(:,2).*in2(:,4).*sin(in2(:,3)).*6.0).*
 (1.0e1./3.0))./(cos(in2(:,3)).^2.*1.5e2-1.87e2)]

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

[Star Strider]

First, change the matlabFunction call to:

ODESYS = matlabFunction(odes, 'Vars',{t [x1 x2 x3 x4] u1, u2});

then use ‘ODESYS’ as:

[T,X] = ode45(@(T,X) ODESYS(T, X, u1, u2), [0 5], x0);

and provide some value for ‘u1’ and ‘u2’ (or rewrite your function to omit them entirely).

Answer 2

[Ajay Pratap Singh]

syms TU(t) TL(t) TP(t)

for i=1:2

% H = [239 527.5];

Ta=[286.15 287.15];

c1=4.10412959805980e-05;

c2=[0.0101769534228370 0.0106028663412749];

c3=3.01522866838092e-06;

c4=1.68909239901093e-05;

c5=[0.00142435898635732 0.00150330661006984];

c6=1.13326159344067e-05;

c7=7.95685343044966e-07;

c8=0.00565004757967287;

c9=[0.000474707322815741 0.00104773268947826];

c10=0.00565004757967287;

TU0(1)=Ta(1);

TL0(1)=Ta(1);

TP0(1)=Ta(1);

eqn1 = (diff(TU,t))+(c1*TU)-(c3*TL)==(c2(i));

% Dr1 = diff(TU,t);

eqn2 = (diff(TL,t))+(c4*TL)-(c6*TP)-(c7*TU)==(c5(i));

% Dr2 = diff(TL,t);

eqn3 = (diff(TP,t))+(c8*TP)-(c10*TL)==(c9(i));

% Dr3 = diff(TP,t);

eqns = [eqn1; eqn2; eqn3]

% V = odeToVectorField(eqns)

cond1 = TU(0) == TU0(i);

cond2 = TL(0) == TL0(i);

cond3 = TP(0) == TP0(i);

% cond1 = [TU(0) == TU0(i), Dr1(0) == TU0(i)];

% cond2 = [TL(0) == TL0(i), Dr2(0) == TL0(i)];

% cond3 = [TP(0) == TP0(i), Dr3(0) == TL0(i)];

conds = [cond1; cond2; cond3]

% [TU(t),TL(t), TP(t)] = dsolve(eqns,conds)

% [TU(t),TL(t), TP(t)] = dsolve(eqns,[0,500],conds)

S = dsolve(eqns,conds)

TU(t) = S.TU;

TL(t) = S.TL;

TP(t) = S.TP;

y1=subs(TU(t),t,3600)

y2=subs(TL(t),t,3600)

y3=subs(TP(t),t,3600)

TU_final=vpa(y1,6)

TL_final=vpa(y2,6)

TP_final=vpa(y3,6)

TU0(i+1)=TU_final(i)

TL0(i+1)=TL_final(i)

TP0(i+1)=TP_final(i)

end

Error using mupadengine/feval_internal

No differential equations found. Specify differential equations by using symbolic functions.

Error in dsolve>mupadDsolve (line 334)

T = feval_internal(symengine,'symobj::dsolve',sys,x,options);

Error in dsolve (line 203)

sol = mupadDsolve(args, options);

Error in Trailforsolarpond (line 37)

S = dsolve(eqns,conds)