Solver: Warning: Unable to find explicit solution

I have the following code:

clear all;

syms a p s e
assume(a>=0)
assume(p>=0 & p<=1)
assume(s<=0)
assume(e>=0)
assumptions
global beta_0 beta_1 pc mc z d b sigma tau phi beta f psi j

beta_0=0.01
beta_1=-0.2
m_e = 1/(1+exp(-(beta_0+beta_1*e)))
m_eprime = diff(m_e, e)
m_eprime2 = diff(m_eprime, e)

pc=-150
mc=-150
z=0.01
d=-5000
b=0
pi = b+z*d+pc+mc

c_a = -a^1.1
c_aprime = diff(c_a, a)
c_aprime2 = diff(c_aprime, a)

sigma=0.88
g = pi+c_a-d  %g>0
sg = s+g  %s+g<0
v_g = g^sigma
v_gprime = sigma*g^(sigma-1)

tau=0.87
phi=0.6
w_p = exp(-tau*(-log(p))^phi)
w_1p = exp(-tau*(-log(1-p))^phi)

beta=0.5
w_beta = exp(-tau*(-log(beta))^phi)
w_1beta = exp(-tau*(-log(1-beta))^phi)

f=20000
df = -d-f
n_h = (a^0.88-e^0.88)^0.5+1
n_heprime = diff(n_h,e)
n_heprime2 = diff(n_heprime,e)
n_heprimea = diff(n_heprime,a)
n_haprime = diff(n_h,a)
n_haprime2 = diff(n_haprime,a)
n_l = (a^0.88-e^0.88)^0.5
n_laprime = diff(n_l,a)
n_laprime2 = diff(n_laprime,a)
n_laprimee = diff(n_laprime,e)
n_leprime = diff(n_l,e)
n_leprime2 = diff(n_leprime,e)

q_h = n_h*f-f
q_l = n_l*f-f

l=d-pi

psi=2.25
v_l = -psi*(-l)^sigma

v_sqh = -psi*(-(s+q_h))^sigma
v_sqhprime = psi*sigma*(-(s+q_h))^(sigma-1);
v_sql = -psi*(-(s+q_l))^sigma
v_sqlprime = psi*sigma*(-(s+q_l))^(sigma-1);

O_a = m_e*w_beta*(w_1p*v_gprime*c_aprime+w_p*v_sqhprime*n_haprime*f)+m_e*w_1beta*(w_1p*v_gprime*c_aprime+w_p*v_sqlprime*n_laprime*f)==0
a_sol = solve(O_a,a, 'Real', true, 'IgnoreAnalyticConstraints', true)
simplify(a_sol)

I need "0_a" solved for "a" so that I can get rid of "a" in following specifications using "a_sol" equation instead. The problem is that the solve function looks unable to handle the input. I've tried all of the recomendations depicted in the help center, but the problem persists. Unfortunately, I cannot simplify or specify "0_a" any more because I need it to be quite flexible so that I can try different values for my globals and different specifications of the functions described above. Can anyone help me?

Answer 1

[Divyank]

The issue you are facing is most likely due to the complexity of the equation you are trying to solve. The equation O_a is a nonlinear equation in terms of 'a', which means that it cannot be solved algebraically in a closed form.

However, you can try using numerical methods to find an approximate solution to the equation. One way to do this is to use the 'vpasolve' function in MATLAB, which uses a numerical solver to find an approximate solution to the equation.

Here is an example of how you can use 'vpasolve' to solve the equation O_a for 'a':

% Define the equation to be solved
eqn = m_e*w_beta*(w_1p*v_gprime*c_aprime+w_p*v_sqhprime*n_haprime*f)+m_e*w_1beta*(w_1p*v_gprime*c_aprime+w_p*v_sqlprime*n_laprime*f)==0;

% Solve the equation for 'a' using vpasolve
a_sol = vpasolve(eqn, a, [0, 100], 'RandomStart', true);

% Simplify the solution
a_sol = simplify(a_sol);

Note that the 'vpasolve' function may not always find a solution, or it may find a solution that is not accurate enough for your needs. In this case, you may need to try different starting points or adjust the solver options to get a better result.

Once you have found a solution for 'a', you can substitute it into the other expressions in your code using the 'subs' function. For example, to substitute the solution for 'a' into the expression for 'n_h', you can do:

n_h = subs(n_h, a, a_sol);

I hope this helps!