I have four interconnected PDEs (first order), initial and boundary conditions. Functions are dependent on two variables "mu" (spatial) and "t" (temporal). The system is non-linear (three of the four equations are linear, three of the four equations are time derivatives). NDSolve reports various bugs, but I can't evaluate what's wrong with my system. So I started to "remove" the equations (the functions that solved these equations were replaced by a constant). If this is inappropriate approach, please let me know. With two linear equations (only time differentiations), I got a solution, reporting this error: "NDSolveValue: The given initial conditions were not consistent with the differential-algebraic equations. NDSolve will attempt to correct the values." My first question is, how may I satisfy the initial condition to be sure it is right?
When I add third equation (linear, differentiation in "mu"), it blows up with the following error messages: "NDSolveValue: The given initial conditions were not consistent with the differential-algebraic equations. NDSolve will attempt to correct the values." "LinearSolve: Matrix SparseArray[Automatic,{123,123},0.,<<1>>] is singular." "NDSolveValue: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions." What is wrong in this case?
I also attached screenshot of error messages in the case of solution of three equations (the first error message appears also in the case of solution of two equations).
Equations:
G[mu_, t_] = 4 \[Pi]*frho0[dmu]*r[mu, t]^2*D[r[mu, t], mu];
w[mu_, t_] = 1 + ep[mu, t]/c^2 + p[mu, t]/(frho0[dmu]*c^2);
w[mu_, t_] = 1;
a[mu_, t_] = 1/w[mu, t];
ep[mu_, t_] = k*frho0[dmu]^(\[Gamma] - 1)/(\[Gamma] - 1);
p[mu_, t_] = (\[Gamma] - 1) ep[mu, t]*frho0[dmu];
equt[mu_,
t_] = -a[mu,
t] (4 \[Pi]*r[mu, t]^2*G[mu, t]/w[mu, t]*D[p[mu, t], mu] + (
mu*gr)/frho0[dmu]^2 +
0*(4 \[Pi]*gr)/c^2*p[mu, t]*
r[mu, t]);(*Eq1, rho[mu,t]->frho0[dmu], m[mu,t]\[Rule]mu*)
equt[mu_,
t_] = -a[mu,
t] (4 \[Pi]*r[mu, t]^2*G[mu, t]/w[mu, t]*D[p[mu, t], mu] + (
m[mu, t]*gr)/frho0[dmu]^2 +
0*(4 \[Pi]*gr)/c^2*p[mu, t]*
r[mu, t]);(*Eq1, rho[mu,t]->frho0[dmu]*)
eqrt[mu_, t_] = a[mu, t]*u[mu, t];(*Eq2*)
eqmm[mu_, t_] =
4 \[Pi]*frho0[dmu]*(1 + ep[mu, t]/c^2)*
r[mu, t]^2 D[r[mu, t], mu];(*Eq3, rho[mu,t]->frho0[dmu]*)
eqrhort[mu_, t_] = -a[mu, t]*frho0[dmu]*r[mu, t]^2 D[u[mu, t], mu]/
D[r[mu, t], mu];(*Eq4, rho[mu,t]->frho0[dmu]*)
(*equations and boudnary and initial conditions for 3 equations*)
eqs = {D[u[mu, t], t] == equt[mu, t], D[r[mu, t], t] == eqrt[mu, t],
D[m[mu, t], mu] == eqmm[mu, t]};
bcon = {DirichletCondition[u[mu, t] == 0., mu == dmu],
DirichletCondition[r[mu, t] == 0., mu == dmu],
DirichletCondition[m[mu, t] == dmu, mu == dmu]};
incon = {u[mu, 0] == 0., r[mu, 0] == r0[mu],
m[mu, 0] == fm0[mu]}; {dmu, mumax}
(*solution*)
{fu, fr, fm} =
NDSolveValue[{eqs, incon, bcon}, {u, r, m}, {mu, dmu, mumax}, {t, 0,
0.1}]
Attachments:
