Short Answer
Set
Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}
inside NDSolve will resolve the problem. It's not necessary to set "ScaleFactor" to 1, it just needs to be a not-that-small positive number.
Long Answer
The answer for this problem is hidden in this obscure tutorial. I'll try my best to retell it in an easier-to-understand way.
Let's consider the following simpler initial-boundary value problem (IBVP) of heat equation that suffers from the same issue:
$$\frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2}$$ $$u(0,x)=x(1-x)$$ $$u(t,0)=0,\ \frac{\partial u}{\partial x}\bigg|_{x=1}=0$$ $$t>0,\ 0\leq x\leq 1$$
Clearly $u(0,x)=x(1-x)$ and $\frac{\partial u}{\partial x}\bigg|_{x=1}=0$ is inconsistent. When you solve it with NDSolve / NDSolveValue, ibcinc warning will be spit out:
tend = 1; xl = 0; xr = 1; With[{u = u[t, x]}, eq = D[u, t] == D[u, x, x]; ic = u == (x - xl) (xr - x) /. t -> 0; bc = {u == 0 /. x -> xl, D[u, x] == 0 /. x -> xr};] sol = NDSolveValue[{eq, ic, bc}, u, {t, 0, tend}, {x, xl, xr}]
NDSolveValue::ibcinc
and further check shows the boundary condition (b.c.) $\frac{\partial u}{\partial x}\bigg|_{x=1}=0$ isn't satisfied at all:
Plot[D[sol[t, x], x] /. x -> xr // Evaluate, {t, 0, tend}, PlotRange -> All]
Why does this happen?
The best way for explaination is re-implementing the method used by NDSolve in this case i.e. method of lines, all by ourselves.
As mentioned in the document, method of lines is a numeric method that discretizes the partial diffential equation (PDE) in all but one dimension and then integrating the semi-discrete problem as a system of ordinary differential equations (ODEs) or differential algebraic equations (DAEs). Here I discretize the PDE with 2nd order centered difference formula:
$$f'' (x_i)\simeq\frac{f (x_{i}-h)-2 f (x_i)+f (x_{i}+h)}{h^2}$$
Clear@dx formula = eq /. {D[u[t, x], t] -> u[x]'[t], D[u[t, x], x, x] -> (u[x - dx][t] - 2 u[x][t] + u[x + dx][t])/dx^2} points = 5; dx = (xr - xl)/(points - 1); ode = Table[formula, {x, xl + dx, xr - dx, dx}]
{u[1/4]'[t] == 16 (u[0][t] - 2 u[1/4][t] + u[1/2][t]), u[1/2]'[t] == 16 (u[1/4][t] - 2 u[1/2][t] + u[3/4][t]), u[3/4]'[t] == 16 (u[1/2][t] - 2 u[3/4][t] + u[1][t])}
I've chosen a very coarse grid for better illustration. The initial condition (i.c.) should also be discretized:
odeic = Table[ic /. u[t_, x_] :> u[x][t] // Evaluate, {x, xl, xr, dx}]
{u[0][0] == 0, u[1/4][0] == 3/16, u[1/2][0] == 1/4, u[3/4][0] == 3/16, u[1][0] == 0}
We still need to deal with b.c.s. The Dirichlet b.c. doesn't need discretization:
bcnew1 = bc[[1]] /. u[t, x_] :> u[x][t]
u[0][t] == 0
The Neumann b.c. contains the derivative of $x$ so we need to discretize it with the one-sided difference formula:
$$f' (x_n)\simeq \frac{f (x_{n}-2h)-4 f (x_{n}-h)+3 f (x_n)}{2 h}$$
bcnew2 = bc[[2]] /. D[u[t, x_], x_] :> (u[x - 2 dx][t] - 4 u[x - dx][t] + 3 u[x][t])/(2 dx)
2 (u[1/2][t] - 4 u[3/4][t] + 3 u[1][t]) == 0
"OK, 5 unknowns, 5 equations, we can now solve the system with any ODE solver! Just as NDSolve does!" Sadly you were wrong if you thought this statement is correct, because:
Though {ode, odeic, bcnew1, bcnew2} is already a solvable system, it's not a set of ODEs, but DAEs. Notice here ODE refers to explicit ODE i.e. the coefficient of derivative term can't be $0$. Clearly, bcnew1 and bcnew2 don't explicitly contain the derivative of $t$.
Though NDSolve is able to handle the DAE system {ode, odeic, bcnew1, bcnew2} directly, it doesn't solve the PDE in this way by default. Instead, it'll try to transform the DAE system to an explicit ODE system, probably because its ODE solver is generally stronger than the DAE solver (at least now).
Remark
To be precise, internally NDSolve doesn't "transform the DAE system to an explicit ODE system", it directly discretize the PDE to an explicit ODE system that can be viewed as a transformation of {ode, odeic, bcnew1, bcnew2}, i.e. {ode, odeic, bcnew1, bcnew2} isn't really generated, but this truth doesn't influence the discussion in this answer.
So, how does NDSolve transform the DAE system to ODE system? "That's simple! Just eliminate some of the variables in ode and odeic with bcnew1 and bcnew2! " Yeah this is a possible method, but not the one implemented in NDSolve. NDSolve has chosen a method that may be rather unusual at first glance. It mixes the original b.c. with its 1st order derivative respects to $t$. For our specific problem, the b.c. becomes:
odebc1 = D[#, t] + scalefactor1 # & /@ bcnew1
scalefactor1 u[0][t] + u[0]'[t] == 0
odebc2 = D[#, t] + scalefactor2 # & /@ bcnew2
2 scalefactor2 (u[1/2][t] - 4 u[3/4][t] + 3 u[1][t]) + 2 (u[1/2]'[t] - 4 u[3/4]'[t] + 3 u[1]'[t]) == 0
Where scalefactor1 and scalefactor2 are properly chosen coefficients.
It's not hard to notice this approach is systematic and easy to implement, and I guess that's the reason why NDSolve chooses it for transforming algebraic equations to ODEs. Nevertheless, this method has its disadvantage. The generated b.c. is equivalent to the original b.c., only if the original b.c. is continuous and i.c. is consistent with the b.c..
Let's use odebc1 as an example. In our case, bc[[1]] is continuous, and it's consistent with ic, so it can be easily rebuilt from odebc1:
DSolve[{odebc1, odeic[[1]]}, u[0][t], t]
{{u[0][t] -> 0}}
However, if the i.c. is something that isn't consistent with bc[[1]], for example u[0][0] == 1, the b.c. rebuilt from odebc1 will become:
DSolve[{odebc1, u[0][0] == 1}, u[0][t], t]
{{u[0][t] -> E^(-scalefactor1 t)}}
It's no longer equivalent to bc[[1]], but when scalefactor1 is a large positive number, this b.c. will converge to the original one.
Now here comes the key point. As stated in the document:
With the default "ScaleFactor" value of Automatic, a scaling factor of 1 is used for Dirichlet boundary conditions and a scaling factor of 0 is used otherwise.
i.e. scalefactor2 will be set to 0. Guess what b.c. will be rebuilt in this way?:
With[{scalefactor2 = 0}, DSolve[{D[#, t] + scalefactor2 # & /@ bc[[2]], D[ic, x] /. x -> xr}, D[u[t, x], x] /. x -> xr, t]]
{{Derivative[0, 1][u][t, 1] -> -1}}
It's a completely different b.c..
Back to the problem mentioned in the question, we can analyse its b.c. with the same method as above:
bcInQuestion = (D[u[t, x], x] /. x -> 1) == c1 (u[t, 1] - c2); icInQuestion = u[0, x] == c3; With[{sf = 0}, DSolve[{D[#, t] + sf # & /@ bcInQuestion, D[icInQuestion, x] /. x -> 1}, D[u[t, x], x] /. x -> 1, t]] /. u[0, _] :> c3
{{Derivative[0, 1][u][t, 1] -> -c1 c3 + c1 u[t, 1]}}
We see c1 is still in the rebuilt b.c., while c2 is completely killed by the derivation, that's why OP found "editing c1 does affect the graph, while editing c2 does not".
OK, then why does NDSolve choose such a strange setting for scaling factor? The document explains as follows:
There are two reasons that the scaling factor to multiply the original boundary condition is zero for boundary conditions with spatial derivatives. First, imposing the condition for the discretized equation is only a spatial approximation, so it does not always make sense to enforce it as exactly as possible for all time. Second, particularly with higher-order spatial derivatives, the large coefficients from one-sided finite differencing can be a potential source of instability when the condition is included. …
but personally I think this design is just too lazy: why not make NDSolve choose a non-zero scaling factor at least when ibcinc warning pops up and the order of spatial derivatives isn't too high (this is usually the case, the differential order of most of PDEs in practise is no higher than 2 )?
Anyway, now we know how to fix the issue. Just choose a positive scaling factor:
Clear[c2] c1 = -10; (*c2=10;*) c3 = 20; q[t, x] = 100000; heat = ParametricNDSolveValue[{1591920 D[u[t, x], t] == .87 D[u[t, x], x, x] + q[t, x], (D[u[t, x], x] /. x -> 0) == 0, (D[u[t, x], x] /. x -> 1) == c1 (u[t, 1] - c2), u[0, x] == c3}, u, {t, 0, 600}, {x, 0, 1}, c2, Method -> {"MethodOfLines", "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]; Plot[heat[#][300, x] & /@ Range[10, 50, 10] // Evaluate, {x, 0.9, 1}, PlotRange -> All, AxesLabel -> {"x(m)", "T"}]
Now c2 influences the solution.
Young has also solved the problem with the new-in-v10 "FiniteElement" method, I guess it's probably because b.c.s are imposed in a completely different way when "FiniteElement" method is chosen, but I'd like not to talk too much about it given I'm still in v9 and haven't looked into "FiniteElement".
