2
$\begingroup$

I'm trying to solve a system of differential equations that looks like this:

Q = 1;
z[t_] := Sqrt[1/3 (y[t]^2/2 + V[x[t]] + r[t])];
z0 := Sqrt[1/3 (y0^2/2 + x0^(2 n)/(2 n) + r0)];
k = 100*z0;

a[t_] := x[t]^2;
b[t_] := 4*r[t]*(k/g[t]) + y[t];
c[t_] := x[t] + y[t]*g[t];
d[t_] := 1 + 4*r[t];

A[t_] := {{a[t], b[t]},{c[t], d[t]}};
B[t_] := {{0}, {-Sqrt[x[t]^alpha]}};
J[t_] := {{j1[t], j2[t]}, {j3[t], j24[t]}};

sys = {x'[t] == y[t]/z[t],
   y'[t] == -3 (1 + Q) y[t] - x[t]/z[t],
   r'[t] == -4 r[t] + 3 Q*y[t]^2,
   g'[t] == g[t],
   J'[t] == -A[t].J[t] - J[t].Transpose[A[t]] + B[t].Transpose[B[t]]};

ICs = {x[0] == 15, y[0] == -0.8, r[0] == 0.005, g[0] == 1, 
   J[0] == {{0, 1}, {1, 0}}};

soln = NDSolve[{sys, ICs}, {x[t], y[t], r[t], J[t], g[t]}, {t, 0, 10},
   Method -> {"StiffnessSwitching"}]

The difference between this and my real code is that, A and J must be 4x4 matrices. I have no idea how to resolve this. Mathematica returns the error:

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`.

But I can't figure out how to write such a problem for 2x2 or 4x4 matrices.

-------- Update ---------

In small versions the suggestion worked great, however for my problem:

    ClearAll["Global`*"]
    Mp = 1; 
    gstar = 228.75; 
    Cr := gstar*(Pi^2/30); 
    lambda = 1; 
    V0 := 1; 
    n = 1; 
    c = 0.012;

   V[x_] := x^(2*n)/(2*n); 
   dV[x_] := D[V[x], x]; 
   ddV[x_] := D[dV[x], x]; 
   T[t_] := (r[t]/Cr)^(1/4); 
  
   Q = 1; 
   z[t_] := Sqrt[(1/3)*(y[t]^2/2 + V[x[t]] + r[t])]; 
   z0 := Sqrt[(1/3)*(y0^2/2 + x0^(2*n)/(2*n) + r0)]; 
   k = 100*z0; 
   beta = 1; 
   alpha = 0; 
   gamma[t_] := c*x[t]^alpha*T[t]^beta; 
   gammaT[t_] := c*x[t]^alpha*beta*T[t]^(beta - 1); 
   gammaPhi[t_] := c*alpha*x[t]^(alpha - 1)*T[t]^beta; 


   fpsi[t_] := 1 + k^2/(3*1^2*z[t]^2) - y[t]^2/6; 
   frho[t_] := 1/(6*z[t]^2); 
   fdphi[t_] := y[t]/6; 
   fphi[t_] := x[t]/(6*z[t]^2); 
   gpsi[t_] := (c*x[t]^alpha*T[t]^beta)*z[t]*y[t]^2 - 
  (k^2/(3*1^2))* ((2*k^2)/(a[t]^2*z[t]^2)) - y[t]^2; 

   grho[t_] := 4 - (((c*x[t]^alpha*beta*T[t]^(beta - 
   1))*z[t]*T[t])/(4*r[t]))*y[t]^2 - k^2/(3*1^2*z[t]^2); 

   gphi[t_] := (-(k^2/(3*a[t]^2*z[t]^2)))*(3*z[t]^2*y[t] + (dV[x] 
   /.{x -> x[t]})) - z[t]*(c*alpha*x[t]^(alpha - 
   1)*T[t]^beta)*y[t]^2;


   gdphi[t_] := (-(k^2/(3*1^2) + 2* 
  (c*x[t]^alpha*T[t]^beta)*z[t]))*y[t];

  hpsi[t_] := (2*x[t]^3)/z[t]^2 + 
  ((c*x[t]^alpha*T[t]^beta)/z[t])*y[t]; 

  hrho[t_] := ((T[t]*(c*x[t]^alpha*beta*T[t]^(beta - 
   1)))/(4*r[t]*z[t]))*y[t]; 

  hdphi[t_] := 3 + (c*x[t]^alpha*T[t]^beta)/z[t] + (1/z[t])*((1/2)* 
   (x[t]^2/2/x[t]^2)); 

   hphi[t_] := k^2/(z[t]^2*1^2) + (ddV[x] /. {x -> x[t]})/z[t]^2 + 
  ((c*alpha*x[t]^(alpha - 1)*T[t]^beta)/z[t])*y[t];

Matrices are formed with these components:

  A[t_] := {{fpsi[t], frho[t], fdphi[t], fphi[t]}, 
  {gpsi[t] + 4*r[t]*fpsi[t], grho[t] + 4*r[t]*frho[t],
   gdphi[t] + 4*r[t]*fdphi[t], gphi[t] + 4*r[t]*fphi[t]}, 
   {hpsi[t] + 4*y[t]*fpsi[t], hrho[t] + 4*y[t]*frho[t], hdphi[t] + 
    4*y[t]*fdphi[t], hphi[t] + 4*y[t]*fphi[t]}, {0, 0, -1, 0}}; 


  B[t_] := {{0}, {(-Sqrt[2*(c*x[t]^alpha*T[t]^beta)*z[t]* 
  (T[t]/1^3)])*y[t]}, {Sqrt[2*(c*x[t]^alpha*T[t]^beta)* 
  (T[t]/(a[t]*z[t])^3)]}, {0}}; 

   J[t_] := {{j11[t], j12[t], j13[t], j14[t]}, {j21[t], j22[t], 
   j23[t], j24[t]}, {j31[t], j32[t], j33[t], j34[t]}, {j41[t], 
   j42[t], j43[t], j44[t]}}; 

   sys = {x'[t] == y[t]/z[t],
    y'[t] == -3 (1 + Q) y[t] - x[t]/z[t],
    r'[t] == -4 r[t] + 3 Q*y[t]^2,
    a'[t] == a[t],
    J'[t] == -A[t].J[t] - J[t].Transpose[A[t]] + 
    B[t].Transpose[B[t]]};

    ICs = {x[0] == 15,
     y[0] == -0.8,
     r[0] == 0.005, 
     a[0] == 1, 
    j11[0] = 0, j12[0] = 0, j13[0] = 0, j14[0] = 0, j21[0] == 0, 
    j22[0] == 0, j23[0] == 0, j24[0] == 0, j31[0] == 0, 
    j32[0] == 0, j33[0] == 1 + I*100, j34[0] == -1 + I*100, 
    j41[0] == 0, j42[0] == 0, j43[0] == 1 + 100, j44[0] == 1}; 


     Block[{n, r0, x0, y0, alpha, beta}, {n, r0, x0, y0, alpha, 
     beta} = {1, 0.005, 15, -0.8, 0, 1}; 
     soln = NDSolve[{sys, ICs}, Flatten[{x[t], y[t], r[t], J[t], 
     a[t]}], 
     {t, 0, 10}, Method -> {"StiffnessSwitching"}]]

I received the error:

 Equation or list of equations expected instead of 0 in the first \
 argument.

And I don't know how to solve it. It must be something simple but I can't see it.

-------- UPDATE 2 ----------

The problem was the part where I wrote 'V[x[t]]', instead I directly write x[t]^2/2, and in the parts where dV[x] appears I write x[t] etc. ... Now it's settled. Thanks.

$\endgroup$
5
  • $\begingroup$ You seemed to have some misplaced backticks in your code, which I tried to fix. You might want to check my edit.. You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is useful for learning how to format your questions and answers. You may also find the meta Q&A, How to copy code from Mathematica so it looks good on this site, helpful $\endgroup$
    – Michael E2
    Aug 26, 2023 at 13:40
  • $\begingroup$ The following checks your system and shows that some of the functions are parameters are undefined (at least in the posted code): Solve[sys, D[Flatten@{x[t], y[t], r[t], J[t], g[t]}, t]] /. t -> 0 /. Solve[ICs, Flatten@{x[t], y[t], r[t], J[t], g[t]} /. t -> 0] $\endgroup$
    – Michael E2
    Aug 26, 2023 at 13:44
  • $\begingroup$ Oh thanks, it was my first question, I'll look at this link, it will be very helpful thanks, $\endgroup$ Aug 26, 2023 at 13:45
  • $\begingroup$ Yes it's true, when writing the example I forgot to add some things like n = 1; z[t_] := Sqrt[1/3 (y[t]^2/2 + x[t] + r[t])]; x0 = 15; y0 = -0.8; r0 = 0.0005;z0 := Sqrt[1/3 (y0^2/2 + x0^(2 n)/(2 n) + r0)]; But, I believe that this type of system can be solved by writing: NDSolve[sys, D[Flatten@{x[t], y[t], r[t], J[t], g[t]}, t] ] /. t -> 0 /. Solve[ICs, Flatten@{x[t], y[t], r[t], J[t], g[t]} /. t -> 0] $\endgroup$ Aug 26, 2023 at 13:51
  • $\begingroup$ If those particular value cause a problem for NDSolve, then you should edit the question to include them. Otherwise, I think my answer will show the correct approach. $\endgroup$
    – Michael E2
    Aug 26, 2023 at 13:53

1 Answer 1

4
$\begingroup$

If I fill in arbitrarily chosen values for the missing parameters and functions, NDSolve computes a solution. Note the use of Flatten, which is not necessary for success, but it produces a more useful output. (Personally, I'd leave off the [t] and solve for {x, y,...}, but I like that form of solution better than the x[t] ->... form.)

Block[{n, r0, x0, y0, alpha, V},
 {n, r0, x0, y0, alpha, V} = {2, 1, 1, 1, 1, # &};
 soln = NDSolve[{sys, ICs}, 
   Flatten@{x[t], y[t], r[t], J[t], g[t]}, {t, 0, 10}, 
   Method -> {"StiffnessSwitching"}]
 ]

Mathematica graphics

$\endgroup$
1
  • $\begingroup$ Oh thank you very much, I believe this will help me a lot. I'll test it in my real code and any updates I'll comment here. $\endgroup$ Aug 26, 2023 at 13:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.