How do I solve an equation with a vector term?

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

Answer 1

[Basil C.]

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.

[Efaz Ejaz]

Hello, thank you for your answer.

I copypasted your code on to my MATLAB script and I'm afraid it doesn't work. It just prints the error message continuously and doesn't generate a workspace, sort of like an infinite loop.

Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve. 
> In solve (line 304)
  In ROCKETMAN (line 18) 
Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve. 

I was really rooting for your code to work. But thanks anyway.

[Walter Roberson]

You will not be able to get a closed form solution for that, so there is no benefit in using solve() to ask for a closed form algebraic-number solution. Instead just change the solve() to vpasolve()

However... there are multiple solutions. Between b = 0 and approximately 74.2, there are three real-valued solutions, two of which are negative and one positive. With larger b, there is one real-valued solution that is positive. There are an unknown number of complex-valued solutions (at least two) as well. You need to decide what you are expecting.

[Basil C.]

clear all
close all
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
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)= vpasolve(eqn1, t);
end
%%
for i=1:numel(b)
eqn2 = -50000 + h_b(i)+v_b(i)*(t-b(i))-0.5*32.2*(t-b(i))^2;
k=  solve(eqn2, t);
T2(1,i)=k(1);
T2(2,i)=k(2);
end

I hope you had made changes to T2 as well.

[Walter Roberson]

I recommend adding an initial condition to vpasolve() as otherwise you could get any of the three solutions.

[Efaz Ejaz]

I am looking to get two separate arrays T1 and T2, both consisting of only positive real numbers. Thanks.

[Efaz Ejaz]

Sorry for the pain but I just wanted to add that I want to be able to find the difference between the values in T2 and the corresponding values in T1.

So, say when b is 100, I'll get one solution for t in eqn1 (let's call it t1). When b is 100, there'll be a particular value for h_b and v_b and I should get two solutions for t in eqn2(t21 and t22). I want to discard the first solution for t in eqn2 (t21), and do (t22-t1) and put all those differences into another array.

I was hoping to do this for all the values of b.

Answer 2

[Walter Roberson]

%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.

[Efaz Ejaz]

Index exceeds the number of array elements (0).

Error in sym/subsref (line 900)
            R_tilde = builtin('subsref',L_tilde,Idx);

Error in Rocket (line 20)
T2(1,i)=k(1);

This is the error I get.

[Walter Roberson]

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)

[Efaz Ejaz]

Thanks a lot, mate. I really do appreciate the help.