Solving system of equations

2 次查看（过去 30 天）

Nrmn 2020-2-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/507406-solving-system-of-equations

评论： Nrmn 2020-3-20

Hi,

I am trying to solve a system of equations. This system is comprised of 4 first-order differential equations and 4 analytical equations, I have 8 unknown variables. Each equation is dependent on at least 2 different variables. Is there a way to solve such a system of equations? I know of the bvp4c function that I could use for the differential equations because I know the boundary conditions. But in order to solve these, I need to include the analytical equations somehow. Any ideas?

Thanks!

12 个评论
显示 10更早的评论隐藏 10更早的评论

Nrmn 2020-2-27

I LEFT THEM OUT

because the system is quite extensive. I thought there might be a general way to procede. Ok. I'll give it a shot and show the equations. These are equations for a plasma device where you must differentiate between 3 plasma species: neutrals (index n), ions (index i), electrons (index e). Neutrals and ions are sometimes combined to heavy particles (index h)

The unknown parameters are: M_h, V, u_e, u_h, T_h, T_e, n_n, n_e

The system is to be solved in 1 spacial dimension x between x = 0 and x = 0.075e-3.

At x = 0, M_h, V, u_e, u_h, T_h, T_e, n_n and n_e are known. At x = 0.075e-3, M_h = 1 is known. Throughout the domain, M_e = u_e/sqrt(gamma*T_e/m_e) = const is valid, with gamma and m_e as constants.

The first Differential equation is the progression of the heavy particles' Mach number:

dM_h/dx = M_h*(1+delta*M_h^2)/((1+M_h^2)m_e*n_e*u_h*A)*((1+gamma*M_h^2)/2*W_h-gamma*M_h^2*(1+gamma*M_h^2)/u_h^2*X_h+(1+gamma*M_h^2)/(2*h_0h)*Y_h)

where delta, gamma, m_e, A are constants.

X_h = (R_ne-m_h*u_h*n_dot)*A,

R_ne = -n_e*m_e*nu_en*(u_h-u_e),

nu_en = 6.6e-19*((T_e/4-1)/(1+(T_e/4)^1.6))*n_n*sqrt(8*q*T_e/(pi*m_e)),

n_dot = n_e*n_n*f(T_e)

The second differential equation is for the plasma potential V:

dV/dx = m_e*(nu_ie+nu_en)/q*(u_e-u_h)+1/(q*n_e*A)*(n_e*k_B*A*dT_e/dx+k_B*T_e*A*dn_e/dx)

where m_e, q, k_B, A are constants,

nu_ie = 2*9e-12*n_e/T_e^(3/2)*(23-ln(10^-6*n_e/T_e^3)),

nu_en = 6.6e-19*((T_e/4-1)/(1+(T_e/4)^1.6))*n_n*sqrt(8*q*T_e/(pi*m_e))

The third and fourth differential equations are as follows:

d/dx(h_0e*m_e*n_e*u_e*A) = Y_e

d/dx(h_0h*m_h*(n_n+n_e)*u_n*A) = Y_h

with

h_0e = gamma/(gamma-1)*q*T_e/m_e+u_e^2/2

h_0h = gamma/(gamma-1)*k_B*T_h/m_h+u_h^2/2

The first analytical equation is the conservation of current:

I_d = q*n_e*A*(ue-u_h) = const

The second analytical equation is the conservation of mass flow:

m_dot = A*m_h*u_h*(n_e+n_n) = const

The third analytical equation is for the electron temperature T_e:

(D_ins/2/B_01)^2*n_n*sigma_ion*sqrt(8*q*T_e/(pi*m_e))-q/m_h*(k_B/q*T_h+T_e)/(sigma_cex*n_n)*sqrt(m_h/(k_B*T_h)) = 0

where D_ins, B_01, sigma_cex are constants, sigma_ion = f(T_e)

The last equation is for the neutral pressure:

(n_e+n_n)*k_B/T_h = sqrt(0.78*mfr*zeta*T_r*L_st/D_st^4)*133.32

with zeta = f(T_h) and T_r = f(T_h)

I doubt that anyone can help me in detail, but maybe there are some tips on proceeding.

Thanks.

darova 2020-3-2

I mean values ;)

Is there any data? Or it's function?

Nrmn 2020-3-2

I do not have specific values or data for Y. However, I found an expression in literature for

and

. Maybe this is helpful.

with

and

with

I hope this helps. I'm starting to lose track... I really appreciate your help!

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

darova 2020-3-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/507406-solving-system-of-equations#answer_418146

在 MATLAB Online 中打开

main.m

Here is algorithm i choosed:

integrate 3 and 4 equations. Get 2 nonlinear equations
having 6 nonlinear equations (red box) use fsolve to calculate unknown u_e u_h T_h T_e n_n n_e
calculate M_h and V

I places these equations into ode45 function. I got some results. BUt how to know if they are correct?

I choosed Y_h=1 and Y_e=1 (couldn't handle it)

After constructing system of equations (RED BOX) i put there initial conditions

%       u_e     u_h    T_h   T_e         n_n        n_e 
u0 = [18032, 30.0623, 4500, 1.37, 1.6218E+22, 1.4459E+21);
EQNS(u0,1,1)
ans =
   1.0e+03 *
    0.0016
   -0.0008
    0.0098
    0.0000
   -0.0001
   -2.6002
% shouldn't all they be zero?

See attached script

8 个评论
显示 6更早的评论隐藏 6更早的评论

Nrmn 2020-3-7

BC_Solver.m

Ok, I gave it a shot with bvp4c using your code as base. I changed a the equation system a little bit. There are 6 analytic equations:

I basically left

out of these equations because I don't think it stays constant along x. Eq. 5 is the condition that

must stay constant along x. Eq. 6 is the

that I now assume to be constant. The constant

is simple derived from initial conditions.

The differential equation solver bvpfun is now only a function of

and V. There are two main differential equations that should be solved:

For

in Eq. 7, we need another differential equation:

I tried to approximate Eq. 9 in line 142.

Okay, I now encounter a problem with my boundary conditions. Basically, I know, that

and

. The results of the code, however do not really match the boundary conditions. Moreover, the variables

do not seem to change along x, which is wrong. Could you take a look at the code and tell me where my thinking is false. You're a saint by the way. Thank you so much.