MATLAB: Solve/fsolve system of equations with exp, and multiple variables.

expfsolvesolvesystem of equations

Hello, wanted to solve system of equations, but cant solve and fsolve don't give answers. Solve:
clc
clear all
syms t tf x(t) x1(t) x1_(t) x2(t) x2_(t) L1(t) L2(t) L1_(t) L2_(t) u(t) u_(t) m k c f x1x2(t) S(t) eq_u(t) L1L2(t) eq_u_C(t) x_C(t) d_H_tf_C(t);
F1=subs(x2_+c/m*x2+k/m*x1,[c k m],[2 10 1])
H=subs(u^2+L1*x2+L2*(-k/m*x1-c/m*x2+u),[c k m],[2 10 1])
H1=subs(H,[x1_ x2_ L1_ L2_],[diff(x1,t) diff(x2,t) diff(L1,t) diff(L2,t)])
F=subs(F1,[x1 x1_ x2 x2_],[x diff(x,t) diff(x,t) diff(x,t,2)])
S=0
d_H_x1=subs(subs(diff(subs(H,x1(t),f),f),f,x1(t)),[x1 x1_ x2 x2_],[x diff(x,t) diff(x,t) diff(x,t,2)])
d_H_x2=subs(subs(diff(subs(H,x2(t),f),f),f,x1(t)),[x1 x1_ x2 x2_],[x diff(x,t) diff(x,t) diff(x,t,2)])
d_H_u=subs(subs(diff(subs(H,u(t),f),f),f,u(t)),[x1 x1_ x2 x2_],[x diff(x,t) diff(x,t) diff(x,t,2)])
d_H_tf(t)=subs(diff(S,t)+H,[x1 x1_ x2 x2_],[x diff(x,t) diff(x,t) diff(x,t,2)])
eq_u(t)=solve(subs(d_H_u,u,f)==0,f)
L1L2=dsolve([d_H_x1==-diff(L1,t) d_H_x2==-diff(L2,t)])
eq_u_C(t)=subs(eq_u,[L1 L2],[L1L2.L1 L1L2.L2])
x_C(t)=dsolve(F==eq_u_C)
d_H_tf_C(t)=subs(d_H_tf,[u L1 L2],[eq_u_C L1L2.L1 L1L2.L2])
C1234=solve([d_H_tf_C(tf)==0 x_C(0)==3 x_C(tf)==3 subs(diff(x_C,t),t,0)==0 subs(diff(x_C,t),t,tf)==0])
solve: Unable to find explicit solution.
Also trying fsolve:
function F = myfun(x)
F=[((x(1)*cos(3*x(5))*exp(x(5)))/2 - (x(2)*sin(3*x(5))*exp(x(5)))/2)^2 + (x(1)*(cos(3*x(5))*exp(x(5)) + 3*sin(3*x(5))*exp(x(5))) + x(2)*(3*cos(3*x(5))*exp(x(5)) - sin(3*x(5))*exp(x(5))))*diff(x(x(5)), x(5)) - (x(1)*cos(3*x(5))*exp(x(5)) - x(2)*sin(3*x(5))*exp(x(5)))*(10*x(x(5)) + 2*diff(x(x(5)), x(5)) + (x(1)*cos(3*x(5))*exp(x(5)))/2 - (x(2)*sin(3*x(5))*exp(x(5)))/2);x(3) - (3*x(2))/80 - x(1)/80 - 3;x(3)*cos(3*x(5))*exp(-x(5)) - (sin(3*x(5))*exp(x(5))*(10*x(1) + x(1)*cos(6*x(5)) + 3*x(2)*cos(6*x(5)) + 3*x(1)*sin(6*x(5)) - x(2)*sin(6*x(5))))/240 - x(4)*sin(3*x(5))*exp(-x(5)) + (cos(3*x(5))*exp(x(5))*(x(2)*cos(6*x(5)) - 3*x(1)*cos(6*x(5)) - 10*x(2) + x(1)*sin(6*x(5)) + 3*x(2)*sin(6*x(5))))/240 - 3;- x(1)/8 - x(3) - 3*x(4) ;x(4)*sin(3*x(5))*exp(-x(5)) - (sin(3*x(5))*exp(x(5))*(x(2)*cos(6*x(5)) - 3*x(1)*cos(6*x(5)) - 10*x(2) + x(1)*sin(6*x(5)) + 3*x(2)*sin(6*x(5))))/80 - x(3)*cos(3*x(5))*exp(-x(5)) - 3*x(4)*cos(3*x(5))*exp(-x(5)) - 3*x(3)*sin(3*x(5))*exp(-x(5)) - (sin(3*x(5))*exp(x(5))*(10*x(1) + x(1)*cos(6*x(5)) + 3*x(2)*cos(6*x(5)) + 3*x(1)*sin(6*x(5)) - x(2)*sin(6*x(5))))/240 + (cos(3*x(5))*exp(x(5))*(6*x(1)*cos(6*x(5)) + 18*x(2)*cos(6*x(5)) + 18*x(1)*sin(6*x(5)) - 6*x(2)*sin(6*x(5))))/240 + (sin(3*x(5))*exp(x(5))*(6*x(2)*cos(6*x(5)) - 18*x(1)*cos(6*x(5)) + 6*x(1)*sin(6*x(5)) + 18*x(2)*sin(6*x(5))))/240 - (cos(3*x(5))*exp(x(5))*(10*x(1) + x(1)*cos(6*x(5)) + 3*x(2)*cos(6*x(5)) + 3*x(1)*sin(6*x(5)) - x(2)*sin(6*x(5))))/80 + (cos(3*x(5))*exp(x(5))*(x(2)*cos(6*x(5)) - 3*x(1)*cos(6*x(5)) - 10*x(2) + x(1)*sin(6*x(5)) + 3*x(2)*sin(6*x(5))))/240 ]
end
with
x0=[2;2;2;2;2]
options = optimoptions('fsolve','Display','iter');
[x,fval]=fsolve(@(x) myfun(x),x0,options)
but gets errors:
Array indices must be positive integers or logical values.
Error in myfun (line 9)
F=[((x(1)*cos(3*x(5))*exp(x(5)))/2 - (x(2)*sin(3*x(5))*exp(x(5)))/2)^2 + (x(1)*(cos(3*x(5))*exp(x(5)) + 3*sin(3*x(5))*exp(x(5))) + x(2)*(3*cos(3*x(5))*exp(x(5)) - sin(3*x(5))*exp(x(5))))*diff(x(x(5)), x(5)) - (x(1)*cos(3*x(5))*exp(x(5)) - x(2)*sin(3*x(5))*exp(x(5)))*(10*x(x(5)) + 2*diff(x(x(5)), x(5)) + (x(1)*cos(3*x(5))*exp(x(5)))/2 - (x(2)*sin(3*x(5))*exp(x(5)))/2);x(3) - (3*x(2))/80 - x(1)/80 - 3;x(3)*cos(3*x(5))*exp(-x(5)) - (sin(3*x(5))*exp(x(5))*(10*x(1) + x(1)*cos(6*x(5)) + 3*x(2)*cos(6*x(5)) + 3*x(1)*sin(6*x(5)) - x(2)*sin(6*x(5))))/240 - x(4)*sin(3*x(5))*exp(-x(5)) + (cos(3*x(5))*exp(x(5))*(x(2)*cos(6*x(5)) - 3*x(1)*cos(6*x(5)) - 10*x(2) + x(1)*sin(6*x(5)) + 3*x(2)*sin(6*x(5))))/240 - 3;- x(1)/8 - x(3) - 3*x(4) ;x(4)*sin(3*x(5))*exp(-x(5)) - (sin(3*x(5))*exp(x(5))*(x(2)*cos(6*x(5)) - 3*x(1)*cos(6*x(5)) - 10*x(2) + x(1)*sin(6*x(5)) + 3*x(2)*sin(6*x(5))))/80 - x(3)*cos(3*x(5))*exp(-x(5)) - 3*x(4)*cos(3*x(5))*exp(-x(5)) - 3*x(3)*sin(3*x(5))*exp(-x(5)) - (sin(3*x(5))*exp(x(5))*(10*x(1) + x(1)*cos(6*x(5)) + 3*x(2)*cos(6*x(5)) + 3*x(1)*sin(6*x(5)) - x(2)*sin(6*x(5))))/240 + (cos(3*x(5))*exp(x(5))*(6*x(1)*cos(6*x(5)) + 18*x(2)*cos(6*x(5)) + 18*x(1)*sin(6*x(5)) - 6*x(2)*sin(6*x(5))))/240 + (sin(3*x(5))*exp(x(5))*(6*x(2)*cos(6*x(5)) - 18*x(1)*cos(6*x(5)) + 6*x(1)*sin(6*x(5)) + 18*x(2)*sin(6*x(5))))/240 - (cos(3*x(5))*exp(x(5))*(10*x(1) + x(1)*cos(6*x(5)) + 3*x(2)*cos(6*x(5)) + 3*x(1)*sin(6*x(5)) - x(2)*sin(6*x(5))))/80 + (cos(3*x(5))*exp(x(5))*(x(2)*cos(6*x(5)) - 3*x(1)*cos(6*x(5)) - 10*x(2) + x(1)*sin(6*x(5)) + 3*x(2)*sin(6*x(5))))/240 ]
Error in cw7>@(x)myfun(x)
Error in finDiffEvalAndChkErr
Error in finitedifferences
Error in computeFinDiffGradAndJac
Error in levenbergMarquardt (line 100)
    [JAC,~,~,numEvals,evalOK] = computeFinDiffGradAndJac(XOUT,funfcn,confcn,costFun, ...
Error in fsolve (line 417)
        levenbergMarquardt(funfcn,x,verbosity,options,defaultopt,f,JAC,caller, ...
Error in cw7 (line 28)
[x,fval]=fsolve(@(x) myfun(x),x0,options)

Best Answer

Your d_H_tf_C(tf)==0 turns out to be a differential equation involving x(tf) . You can solve the remaining 4 equations to derive C1, C2, C3, C4, and substitute that in to d_H_tf_C(tf)==0, leaving you with a single differential equation in terms of tf and x(tf) .
When solve() sees a differential equation, it automatically invokes dsolve() on your behalf, but you could make the steps more clear by doing that solve for C1, C2, C3, C4, and trying the dsolve() on what remains.
Unfortunately what remains is too complicated for MATLAB to be able to solve.
When I pass it over to Maple, Maple generates a solution in terms of two indefinite integrals and an unresolved boundary condition. The two indefinite integrals do not appear to have closed form solutions. 
EQN = -19440*exp(tf)*(-8/81*exp(3*tf)*(10*x(tf)+diff(x(tf),tf)-30)*cos(3*tf)^4+8/27*( ...
sin(3*tf)*exp(3*tf)+3/2*exp(2*tf)-3/2*exp(4*tf))*diff(x(tf),tf)*cos(3*tf)^3+((( ...
-80/9+40/27*x(tf)+4/27*diff(x(tf),tf))*exp(2*tf)-4/9*(-20+10/3*x(tf)+diff(x(tf) ...
,tf))*exp(4*tf))*sin(3*tf)+(520/81*x(tf)+16/81*diff(x(tf),tf)+200/27)*exp(3*tf) ...
+(-20/9*x(tf)+2/9*diff(x(tf),tf)-20/3)*exp(5*tf)-2/9*exp(tf)*(10*x(tf)+diff(x( ...
tf),tf)+30))*cos(3*tf)^2+2/3*diff(x(tf),tf)*((exp(tf)-22/9*exp(3*tf)+exp(5*tf)) ...
*sin(3*tf)-31/6*exp(2*tf)+31/6*exp(4*tf)-3/2*exp(6*tf)+3/2)*cos(3*tf)+((-49/27* ...
diff(x(tf),tf)-310/27*x(tf)+80/9)*exp(2*tf)+(310/27*x(tf)+25/9*diff(x(tf),tf)- ...
80/9)*exp(4*tf)+(-10/3*x(tf)-diff(x(tf),tf))*exp(6*tf)+10/3*x(tf)+1/3*diff(x(tf ...
),tf))*sin(3*tf)+(-280/27+235/81*diff(x(tf),tf)-440/81*x(tf))*exp(3*tf)+(20/3+ ...
20/9*x(tf)-29/9*diff(x(tf),tf))*exp(5*tf)+exp(7*tf)*diff(x(tf),tf)-7/9*exp(tf)* ...
(-60/7-20/7*x(tf)+diff(x(tf),tf)))/(4*exp(2*tf)*cos(3*tf)^2-22*exp(2*tf)+9*exp( ...
4*tf)+9)^2 == 0
and Maple's dsolve() says
x(tf) = (300*int(-2/27*(cos(3*tf)+1)*(cos(3*tf)-1)*(-1/3*exp(2*tf)*cos(3*tf)^2+ ...
(exp(tf)-exp(3*tf))*sin(3*tf)-7/6*exp(2*tf)+3/4*exp(4*tf)+3/4)*exp((2+3*i)*tf)/ ...
(-4/27*cos(3*tf)^4*exp(3*tf)+(4/9*sin(3*tf)*exp(3*tf)+2/3*exp(2*tf)-2/3*exp(4* ...
tf))*cos(3*tf)^3+((2/9*exp(2*tf)-2/3*exp(4*tf))*sin(3*tf)-1/3*exp(tf)+8/27*exp( ...
3*tf)+1/3*exp(5*tf))*cos(3*tf)^2+((exp(tf)-22/9*exp(3*tf)+exp(5*tf))*sin(3*tf)- ...
31/6*exp(2*tf)+31/6*exp(4*tf)-3/2*exp(6*tf)+3/2)*cos(3*tf)+(-49/18*exp(2*tf)+25 ...
/6*exp(4*tf)-3/2*exp(6*tf)+1/2)*sin(3*tf)-7/6*exp(tf)+235/54*exp(3*tf)-29/6*exp( ...
5*tf)+3/2*exp(7*tf))/((3+I)*exp((1+3*i)*tf)-4/3-i-5/3*exp(6*i*tf)),tf) + CONSTANT)* ...
exp(-10*int((2*cos(3*tf)^2*exp(tf)+(3*exp(2*tf)-3)*sin(3*tf)-2*exp(tf))/(2*cos( ...
3*tf)^2*exp(tf)+(-6*sin(3*tf)*exp(tf)+9*exp(2*tf)-9)*cos(3*tf)+(9*exp(2*tf)-3)* ...
sin(3*tf)+7*exp(tf)-9*exp(3*tf)),tf))
vpasolve() cannot deal with the system because it vpasolve() cannot deal with differential equations.
If you had one more boundary condition then you could do numeric simultations using ode45() or related (looks like it might be stiff, so perhaps ode23s). In that regard, see odeFunction()
Related Question