Блог

Я совсем новичок в Джулии, и я рассматриваю следующую проблему. Я хотел бы решить (возможно, жесткую) систему ODE, которая описывает релаксацию потока за ударной волной в соответствии с подходом от состояния к состоянию, что означает, что каждый колебательный уровень молекулярных частиц рассматривается как псевдовидовой с его уравнением непрерывности. Здесь я рассматриваю бинарную смесь N2 / N (фактически концентрация N = 0).

Я разделил код Julia на несколько файлов .jl. В основном я вызываю решатель ODE следующим образом:

 prob = ODEProblem(rpart!,Y0_bar,xspan, 1.)                                                                                          
sol  = DifferentialEquations.solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, save_everystep=true, progress=true)
 

где Y0_bar и xspan были определены ранее, и в файле rpart.jl я определяю систему:

 function rpart!(du,u,p,t)

  ni_b = zeros(l);

  ni_b[1:l] = u[1:l]; print("ni_b = ", ni_b, "n")
  na_b = u[l 1];      print("na_b = ", na_b, "n")
  v_b  = u[l 2];      print("v_b = ",  v_b,  "n")
  T_b  = u[l 3];      print("T_b = ",  T_b,  "n")

  nm_b = sum(ni_b);   #print("nm_b = ", nm_b, "n")
  Lmax = l-1;         #println("Lmax = ", Lmax, "n")
  temp = T_b*T0;      #print("T = ", temp, "n")

  ef_b = 0.5*D/T0;    #println("ef_b = ", ef_b, "n")

  ei_b = e_i./(k*T0); #println("ei_b = ", ei_b, "n")
  e0_b = e_0/(k*T0);  #println("e0_b = ", e0_b, "n")

  sigma   = 2.;                     #println("sigma = ", sigma, "n")
  Theta_r = Be*h*c/k;               #println("Theta_r = ", Theta_r, "n")
  Z_rot   = temp./(sigma.*Theta_r); #println("Z_rot = ", Z_rot, "n")

  M  = sum(m); #println("M = ", M, "n")
  mb = m/M;    #println("mb = ", mb, "n")

  A = zeros(l 3,l 3)

  for i = 1:l
    A[i,i]   = v_b
    A[i,l 2] = ni_b[i]
  end

  A[l 1,l 1] = v_b
  A[l 1,l 2] = na_b

  for i = 1:l 1
    A[l 2,i] = T_b
  end
  A[l 2,l 2] = M*v0^2/k/T0*(mb[1]*nm_b mb[2]*na_b)*v_b
  A[l 2,l 3] = nm_b na_b

  for i = 1:l
    A[l 3,i] = 2.5*T_b ei_b[i] e0_b
  end
  A[l 3,l 1] = 1.5*T_b ef_b
  A[l 3,l 2] = 1/v_b*(3.5*nm_b*T_b 2.5*na_b*T_b sum((ei_b. e0_b).*ni_b) ef_b*na_b)
  A[l 3,l 3] = 2.5*nm_b 1.5*na_b

  AA = inv(A); println("AA = ", AA, "n", size(AA), "n")

  # Equilibrium constant for DR processes
  Kdr = (m[1]*h^2/(m[2]*m[2]*2*pi*k*temp))^(3/2)*Z_rot*exp.(-e_i/(k*temp))*exp(D/temp); println("Kdr = ", Kdr, "n")

  # Equilibrium constant for VT processes
  Kvt = exp.((e_i[1:end-1]-e_i[2:end])/(k*temp)); println("Kvt = ", Kvt, "n")

  # Dissociation processes
  kd = zeros(2,l)
  kd = kdis(temp) * Delta*n0/v0;
  println("kd = ", kd, "n", size(kd), "n")

  # Recombination processes
  kr = zeros(2,l)
  for iM = 1:2
    kr[iM,:] = kd[iM,:] .* Kdr * n0
  end
  println("kr = ", kr, "n", size(kr), "n")

  RD  = zeros(l)

  for i1 = 1:l

    RD[i1] = nm_b*(na_b*na_b*kr[1,i1]-ni_b[i1]*kd[1,i1])   na_b*(na_b*na_b*kr[2,i1]-ni_b[i1]*kd[2,i1])

  end

  println("RD = ",  RD,  "n", size(RD))

  B      = zeros(l 3)
  for i  = 1:l
    B[i] = RD[i]
  end
  B[l 1] = - 2*sum(RD)

  du     = AA*B

end
 

Проблема в том, что когда я запускаю симуляцию и строю решение, похоже, что ничего не произошло, и все профили равны и плоские. Фактически, решения на каждом временном шаге равны самим себе. Итак, я думаю, что допустил какую-то ошибку в обновлении u и du, но я не могу это исправить.
В версии Matlab я получаю правильную эволюцию.

С уважением, Лоренцо

Вы используете версию для изменения выходных данных, но вместо изменения выходных данных вы создаете массив. du .= AA*B