Solving system of equations

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!

[darova]

WHERE

ARE

EQUATIONS?

[Nrmn]

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]

Do you have hose equtions of paper or LaTeX format? It's hard to read this as code

I have a question:

• Is it possible re-write 3d and 4th equations as? You didn't describe Y_e and Y_h. Are they constants?

% 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
(h_0e*m_e*n_e*u_e*A) = Y_e*x
(h_0h*m_h*(n_n+n_e)*u_n*A) = Y_h*x

[Nrmn]

Okay sure:

The unknown parameters are: . The system is to be solved in 1 spacial dimension x between and . At are known. At , is known. Throughout the domain, is valid, with γ and as constants.

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

where are constants, , , ,

and is specified below.

The second differential equation is for the plasma potential V:

where are constants, and

The third and fourth differential equations are as follows:

with

with

and are energy source terms. to be honest, I am not 100 % sure if these stay constant, but let's assume it for now. As you mentined, these equations could thus be rewritten as you stated.

The first analytical equation is the conservation of current:

The second analytical equation is the conservation of mass flow:

The third analytical equation is for the electron temperature :

where are constants,

The last equation is for the neutral pressure:

with and

I appreciate your help!

[darova]

Thank you for explanations and LaTeX formulas

Can you please attach all data you have? All constants, parameters and so on

[darova]

Below are parameters, function and initial conditions necessary for solving

(deal( 1 ) means assigning 'one' to all variable)

% where delta, gamma, m_e, A are constants.
[delta, gamma] = deal( 1 );

% where m_e, q, k_B, A are constants,
[m_e, q, k_B, A] = deal( 1 );

% where D_ins, B_01, sigma_cex are constants, sigma_ion = f(T_e)
[D_ins, B_01, sigma_cex] = deal( 1 );

[M_h, u_e, u_h, T_h, T_e, n_n, n_e] = deal( 1 );
[nu_en, m_h, Y_h] = deal( 1 );
[dT_e, dn_e] = deal( 1 );
[L_st, D_st] = deal( 1 );
f = @(x) 1;

[darova]

I'm stil waiting. Wanna know something

[Nrmn]

Hi,

here are all the constants needed:

, , , , , , , , , , , , ,

And these are all the functions I got:

where ,

,

The Boundary conditions are as follows:

At :

At :

Thanks!

[darova]

Energy?

[Nrmn]

Correct! Or to be more precise, it's an energy source term. It's the energy gain/loss by the respective particle species per length and time. So the unit is W/m.

[darova]

I mean values ;)

Is there any data? Or it's function?

[Nrmn]

I do not have specific values or data for Y. However, I found an expression in literature for and . Maybe this is helpful.

So

with

with

and

with

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

Answer 1

[darova]

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

[Nrmn]

Thanks alot!

I am a little confused to be honest.

What is this vector u? And I recon the input on the right hand side are some kind of initial conditions? What do they represent? And did you specify the boundary conditions on the right end of the computing area (at x = 0.75 mm)

[darova]

u vector is [integral(Y_h) integral(Y_e) M_h V]

du = [Y_h; Y_e; dM_h; dV];

I first wanted to understand if this works so i didn't write boundary conditions at x=0.75mm

DO you understand how ode45 works? Is it clear for you?

Do you have any data to compare with? I have no idea if results are close to solution

Do you understand how system of equations (EQNS) is written? Did you try to pass some values?

[Nrmn]

To be honest.. I am not sure how your script works. Could you elaborate a little bit? I also never worked with neither ode45 nor fsolve. However, I get how you smmarized the equations to EQNS but I don't know what you mean by "passing values". Where do I do that and what is the matrix of output variables?

Concerning data to compare: I do have some data at the exit (x= 0.75 mm). n_n should be around 0.2-0.3e22, n_e is around 0.9e21, V at (x = 0) is 8 (I think I forgot to mention this boundary condition) and at (x = 0.75 mm) 12-13. T_e at the right boundary is around 2.5. These aren't fixed values but the trends should be visible in the results.

[darova]

Simple example of fsolve usage

clc,clear
    % some simple syste of equations
F = @(x) [2*x(1) - 1
        x(2)/2 - 3]
opt = optimset('display','off');
x0 = fsolve(F,[1 1],opt)
    % if i pass (put inside) roots of equations
    % all should be zero
F(x0)

I re-wrote all equations to be solved as follows. If all parameters are correct EQNS should return zeros

Some usefull tips:

• Use %% for creating of sections. Move caret between them and press Ctrl+Enter to run part of code

• Run code until specific line (create breakpoint F12)

Please see how ode45 works (it's main solver for ordinary diff equations)

Generally bvp4c should be used for this case

[Nrmn]

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:

1)

2)

3)

4)

5)

6)

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:

7)

8)

For in Eq. 7, we need another differential equation:

9)

I tried to approximate Eq. 9 in line 142.

Okay, I now encounter a problem with my boundary conditions. Basically, I know, that at , at and at . 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.

[darova]

Only M_h and V are changing

Very simplified version of first six equations ( 1)-6) )

Imagine these lines in a loop (nothing changes after first interation)

x1 = fsolve(@(x)x-2,x0);      % x1 = 2 now
x0 = x1;                      % x0 = 2 now

Analytical equations don't depend on M_h nor V

That is why RES variable doesn't change. In another words

x1 = fsolve(@(x)F(x),x0);           % no M_h or V. So x1 are always the same
dM_h = func(x1, M_h);

I'd approximate Y_h as:

iY_h1 = (h_0h1*m_h*(n_n1+n_e1)*u_h1*A);
iY_h0 = (h_0h0*m_h*(n_n0+n_e0)*u_h0*A);
Y_h = (iY_h1-iY_h0)/dx;

Mh plot V plot

Doesn't look that bad

[darova]

how is it going

[Nrmn]

So I tried to change the system of equation a bit. I now want to solve 4 differential equations:

I also changed the last equation in the analytical system of equations. It is now dependent in .

However, I stil can't get results...