2 views (last 30 days)

How do I solve an equation that has one vector term with the rest being constant values? If you look at the part after "syms t", I am attempting to solve for "t" but there are 1001 values for m_0 in eqn1 and a 1001 values for every value of b, h(b) and v(b). So, I'm hoping to get a vector T1 containing 1001 values of t for 1001 values of m_0 and the same for T2 for every value of b, h(b) and v(b).

clear, clc

%Write the given values, u, m_e, b, q, and g.

u = 8000; m_e = 1500; g = 32.2;q = 15; t_0 = 0; b = 0:0.1:100;

%Compute m_0, h_b, and v_b.

m_0 = m_e + q.*b;

h_b = ((u.*m_e)./q)*log(m_e./(m_e+q.*b))+u.*b - 0.5.*g.*b.^2;

v_b = u*log(m_0/m_e) - g.*b;

%Part I'm not so sure about.

syms t

eqn1 = 50000 == u./q.*(m_0-q.*t).*log(m_0-q.*t)+u.*(log(m_0)+1).*t-0.5.*g.*t.^2-((m_0.*u)./q).*log(m_0);

T1 = solve(eqn1, t);

eqn2 = 50000 == h_b+v_b.*(t-b)-0.5.*32.2.*(t-b).^2;

T2 = solve(eqn2, t);

T2_desired = double(T2(2));

T = T2_desired - T1

The output I get is as follows:

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

> In solve (line 317)

In ROCKETMAN (line 15)

Index exceeds the number of array elements (0).

Error in sym/subsref (line 900)

R_tilde = builtin('subsref',L_tilde,Idx);

Be Matlabi
on 11 Nov 2019

The line

eqn1 = 50000 == u./q.*(m_0-q.*t).*log(m_0-q.*t)+u.*(log(m_0)+1).*t-0.5.*g.*t.^2-((m_0.*u)./q).*log(m_0);

will define 1001*1001 equations, you can try indexing instead.

syms t

for i=1:numel(b)

eqn1 = -50000 + u/q*(m_0(i)-q*t)*log(m_0(i)-q*t)+u*(log(m_0(i))+1)*t-0.5*g*t^2-((m_0(i)*u)/q)*log(m_0(i));

T1(i)= solve(eqn1, t);

end

The solution for each will be stored in T, you can do the same for eqn2. Hope this is what you are looking for.

Walter Roberson
on 11 Nov 2019

Walter Roberson
on 11 Nov 2019

%Write the given values, u, m_e, b, q, and g.

u = 8000; m_e = 1500; g = 32.2;q = 15; t_0 = 0; b = 0:0.1:100;

%Compute m_0, h_b, and v_b.

m_0 = m_e + q.*b;

h_b = ((u.*m_e)./q)*log(m_e./(m_e+q.*b))+u.*b - 0.5.*g.*b.^2;

v_b = u*log(m_0/m_e) - g.*b;

%Part I'm not so sure about.

syms t

nb = numel(b);

T1 = zeros(1,nb);

t0 = 1;

for i=1:nb

eqn1 = -50000 + u/q*(m_0(i)-q*t)*log(m_0(i)-q*t)+u*(log(m_0(i))+1)*t-0.5*g*t^2-((m_0(i)*u)/q)*log(m_0(i));

T1(i)= vpasolve(eqn1, t, t0);

t0 = T1(i);

end

%%

T2 = zeros(2,nb);

T0 = 1;

for i = 1:nb

eqn2 = -50000 + h_b(i)+v_b(i)*(t-b(i))-0.5*32.2*(t-b(i))^2;

k = vpasolve(eqn2, t, t0);

T2(1,i) = k(1);

T2(2,i) = k(2);

t0 = T2(1,i);

end

The above is not all that fast.

Note: b = 37.3 is the first b value for which the T2() are real-valued. With smaller b, both T2 are complex valued.

Walter Roberson
on 11 Nov 2019

Even if you had deleted the comment about not so sure, then the line you indicate would be on line 22, not line 20. You must not have used the same code that I posted, so it is difficult to know what might have happened. Also, you did not mention which MATLAB release you are using.

Here is a modified version that takes into account the possibility of too few solutions:

%Write the given values, u, m_e, b, q, and g.

u = 8000; m_e = 1500; g = 32.2;q = 15; t_0 = 0; b = 0:0.1:100;

%Compute m_0, h_b, and v_b.

m_0 = m_e + q.*b;

h_b = ((u.*m_e)./q)*log(m_e./(m_e+q.*b))+u.*b - 0.5.*g.*b.^2;

v_b = u*log(m_0/m_e) - g.*b;

syms t

nb = numel(b);

T1 = zeros(1,nb);

t0 = 1;

wb = waitbar(0, 'T1');

for i=1:nb

waitbar(i./nb, wb);

eqn1 = -50000 + u/q*(m_0(i)-q*t)*log(m_0(i)-q*t)+u*(log(m_0(i))+1)*t-0.5*g*t^2-((m_0(i)*u)/q)*log(m_0(i));

T1(i)= vpasolve(eqn1, t, t0);

t0 = T1(i);

end

%%

T2 = zeros(2,nb);

T0 = 1;

waitbar(0, wb, 'T2');

for i = 1:nb

waitbar(i./nb, wb);

eqn2 = -50000 + h_b(i)+v_b(i)*(t-b(i))-0.5*32.2*(t-b(i))^2;

k = vpasolve(eqn2, t, t0);

if length(k) >= 1

T2(1,i) = k(1);

t0 = T2(1,i);

else

T2(1,i) = nan;

end

if length(k) >= 2

T2(2,i) = k(2);

else

T2(2,i) = nan;

end

end

delete(wb)

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.