Unable to find explicit solution

Can someone please help me to figure out what is wrong with this code. When I run it I got the following error

Warning: Unable to find explicit solution. For options, see help.

Error using mupadengine/feval_internal

Invalid equation or initial condition.

Error in dsolve>mupadDsolve (line 334)

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

Error in dsolve (line 203)

sol = mupadDsolve(args, options);

sol_h = dsolve(eq1,eq2,eq3,eq4,eq5,eq6);

% state equation 
syms x1(t) x2(t) x3(t) w1(t) w3(t) p1(t) p2(t) p3(t)
DX1 = -w3*x2;
DX2 = w3*x1 - w1*x3;
DX3 = w1*x2;

% Cost function inside the integral
syms g;
g = pi/sqrt(2);

% Hamiltonian 
syms x1(t) x2(t) x3(t) w1(t) w3(t) p1(t) p2(t) p3(t) H
H = g -(p1*w3*x2) + p2(w3*x1-w1*x3)+(p3*w1*x2);

% costate equations 
Dp1 = -diff(H,x1);
Dp2 = -diff(H,x2);
Dp3 = -diff(H,x3);

% solve for control w
dw1 = diff(H,w1);
sol_w1 = solve(dw1,w1);
dw3 = diff(H,w3);
sol_w3 = solve(dw3,w3);

% substitute w to state eq
DX1 = subs(DX1, {w1 w3}, {sol_w1 sol_w3});
DX2 = subs(DX2, {w1 w3}, {sol_w1 sol_w3});
DX3 = subs(DX3, {w1 w3}, {sol_w1 sol_w3});

% convert symbolic objects to string for using 'dsolve'
eq1 = strcat('DX1=',char(DX1));
eq2 = strcat('DX2=',char(DX2));
eq3 = strcat('DX3=',char(DX3));
eq4 = strcat('Dp1=',char(Dp1));
eq5 = strcat('Dp2=',char(Dp2));
eq6 = strcat('Dp3=',char(Dp3));
sol_h = dsolve(eq1,eq2,eq3,eq4,eq5,eq6);

% Use Boundary condition
syms cond1 cond2 cond3 cond4 cond5 cond6
cond1 = 'x1(0) == 1';
cond2 = 'x2(0) == 0';
cond3 = 'x3(0) == 0';
cond4 = 'x1(2  == 0';
cond5 = 'x2(2) == 0';
cond6 = 'x3(2) == 1';
sol_H = dsolve(eq1,eq2,eq3,eq4,eq5,eq6,cond1,cond2,cond3,cond4,cond5,cond6);

Thank you for your help.

Answer 1

[Walter Roberson]

In R2020a and earlier, you cannot differentiate with respect to a function. The ability to differentiate with respect to a function was added in R2020b. Here is replacement code for the first part:

Q = @(v) sym(v);

% state equation 
syms x1(t) x2(t) x3(t) w1(t) w3(t) p1(t) p2(t) p3(t)
syms X1 X2 X3 W1 W3 P1 P2 P3
DX1 = -w3*x2;
DX2 = w3*x1 - w1*x3;
DX3 = w1*x2;

% Cost function inside the integral
syms g;
g = Q(pi)/sqrt(Q(2));

% Hamiltonian 
syms x1(t) x2(t) x3(t) w1(t) w3(t) p1(t) p2(t) p3(t) H
H = g -(p1*w3*x2) + p2(w3*x1-w1*x3)+(p3*w1*x2);

% costate equations 
Dp1 = subs(-diff( subs(H,x1,X1), X1), X1, x1);
Dp2 = subs(-diff( subs(H,x2,X2), X2), X2, x2);
Dp3 = subs(-diff( subs(H,x3,X3), X3), X3, x3);

% solve for control w
dw1 = subs(diff( subs(H, w1, W1), W1), W1, w1);
sol_w1 = subs(solve( subs(dw1, w1, W1), W1), W1, w1);
dw3 = subs(diff( subs(H,w3,W3), W3), W3, w3);
sol_w3 = subs(solve( subs(dw3, w3, W3), W3), W3, w3);

However, sol_w1 and sol_w3 both come out empty. dw1 is

p3(t)*x2(t) - D(p2)(w3(t)*x1(t) - w1(t)*x3(t))*x3(t) where that D(p2)(expression) indicates that the derivative of p2 is to be evaluated at the expression. As that expression includes w1 and you want to solve with respect to w1, that would be equivalent to asking which function w1 has the property that the derivative of the unknown function p2 evaluated at that expression, and multiplied by x3(t), comes out equal to p3(t)*x2(t) -- not something you should be able to expect solve() by itself to figure out, and not something that you should expect dsolve() to be able to figure out without further information about the other functions.

Your H is

H = g -(p1*w3*x2) + p2(w3*x1-w1*x3)+(p3*w1*x2);

Notice that you have function p2 being invoked at the difference of two pairs of functions being multiplied together. That's a complicated situation to be able to correctly take the derivative of.

Is it possible that your H should be

H = g -(p1*w3*x2) + p2*(w3*x1-w1*x3)+(p3*w1*x2);

?

[Kaydian Quintero]

Thank you so much.

I did made a mistake with the Hamiltonian (H).

I dont understand the function handle that you created. Can you please explain fuurther.

Again thank you so much

[Walter Roberson]

Q = @(v) sym(v);

accepts a number (such as a floating point number) and asks the symbolic engine to find the best representation it can for the number as one of

• a small integer multiple of pi or pi divided by a small integer
• an integer
• square root of a rational
• a rational

In particular you have

g = pi/sqrt(2);

and I wrap that to

g = Q(pi)/sqrt(Q(2));

Q(pi) -> sym(pi) -> the internal representation of the infinitely long transcendental constant Pi

Q(2) -> sym(2) -> symbolic integer 2

sqrt(Q(2)) -> sqrt() of symbolic integer 2 -> 2^(1/2)

giving an overall result of Pi / 2^(1/2) .

When you do not guide the symbolic engine in this way, the the pi/sqrt(2) is calculated in double precision, and when that double precision is later encounted by the symbolic engine, because it is not able to resolve it to any of the first three cases I list, it ends up converting the number to the rational 1250560371546297 / 562949953421312 which is not very informative.