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

5 visualizaciones (últimos 30 días)
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)
  2 comentarios
darova
darova el 23 de Mayo de 2020
Try vpasolve instead of solve
RobaczeQ
RobaczeQ el 24 de Mayo de 2020
Vpasolve is used when solve cant find results.

Iniciar sesión para comentar.

Respuesta aceptada

Walter Roberson
Walter Roberson el 23 de Mayo de 2020
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()
  9 comentarios
Walter Roberson
Walter Roberson el 25 de Mayo de 2020
syms C1 C2
syms t tf nonnegative
F = ((C2*sin(3*t)*exp(t))/2 - (C1*cos(3*t)*exp(t))/2)^2;
Fi = simplify(int(F, t, 0, tf));
t1 = simplify(diff(Fi,C1));
sol1 = simplify(solve(t1,C1));
Fi2 = simplify(subs(Fi, C1, sol1));
sol1 is optimized by the above code, leaving Fi2 in terms of C1 and tf.
If you examine the parts of Fi2 and do some drawing, you will find that Fi2 is strictly increasing as C2 > 0 increases or as tf increases. The increases is not constant slope in tf, but it is still monotonically increasing. When C2 < 0, Fi2 monotonically decreases toward 0 as C2 increases -- so the minima is when C2 is as close to 0 as you wish to permit.
Therefore, to get the minima, choose C2 and tf arbitrarily close to 0; you can then calculate optimal C1 from sol1. You can reach Fi values arbitrarily close to 0 by choosing small enough (absolute value) C2 and tf.
The analysis shows that there is no magic combination of C2 and tf that are not arbitrarily close to 0 that minimizes Fi. There are just the trivial solutions plus values arbitrarily close to them that do "good enough".
RobaczeQ
RobaczeQ el 26 de Mayo de 2020
Thank You very much for you time and answers. That last two sentences summarise my question.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by