## Issues with no flux boundary conditions in NDSolve

The no flux boundary conditions are making trouble. How to handle them? The problem I’m trying to solve is as follows: PDE: 0.036*D[C[x,y],x]==2.4*10^-10*(D[C[x,y],{x,2}]+D[C[x,y],{y,2}]) Boundary conditions: 1. C[0,y]==Piecewise[{{1,H/21)==0 3. (D[C[x,y],y]/.y->0)==0 4. (D[C[x,y],y]/.y->0.005)==0

I wrote the following: sol2=NDSolve[{0.036*D[C[x,y],x]==2.4*10^-10*(D[C[x,y],{x,2}]+D[C[x,y],{y,2}]),C[0,y]==Piecewise[{{1,H/21)==0,(D[C[x,y],y]/.y->0)==0,(D[C[x,y],y]/.y->0.005)==0},C,{x,0,1},{y,0,0.005}]

The error message shows-The dependent variable in the boundary condition Dirichlet condition need to be linear. But I noticed that a similar no flux boundary condition works just fine (for a different problem though). sol3=NDSolve[{D[u[t, x], t] == D[u[t, x], x, x],u[0, x] == 0,u[t, 0] == Sin[t],(D[u[t, x], x] /. x -> 0)== 0},u, {t, 0, 10}, {x, 0, 5}]

What is the issue in the previous case and how to resolve it?

## Solve boundary value problem with NDSolve. How to print out approximations to a solution?

I solve particular boundary-value-problem for ODE with NDSolve “Shooting” method. The case is that solution is attained very slow, that seems to mean that boundary-value-problem which is supposed to be well-defined enough, in fact, is ill at some place. So i try figure out. First step is to see concrete values of a produced approximations to a solution. What language constructs should i use for that?

Simple example. Suppose we consider particle motion in vertical plane. We throw our particle from initial point with coordinates {x0,y0} := {0,0} and initial trajectory angle 45 degrees. And try to achieve point with coordinates {x1,y1} := {1,0} by varying particle initial velocity. We don’t know two things here: initial velocity and a duration of a motion. Here is how that toy problem can be presented and solved in mathematica:

gravity = 10; bvpsol = NDSolve[     {      {       (* ODE system (5th order) *)       {x''[u] / c[u]^2 == 0,        y''[u] / c[u]^2 == -gravity,        c'[u] == 0},       (* boundary conditions (5 items) *)       {x[0] == y[0] == 0,        x[1] == 1,        y[1] == 0,        x'[0] == y'[0]}       }      }, {x[u], y[u], c[u]}, u,      Method -> {"Shooting",        "StartingInitialConditions" -> {x[0] == y[0] == 0,          x'[0] == y'[0] == 1, c[0] == 1}}] // Flatten;  {dxdu, dydu} = D[{x[u], y[u]} /. bvpsol, u]; {vx0, vy0} = ({dxdu, dydu} / c[u] /. bvpsol) /. {u -> 0}; duration = c[u] /. bvpsol /. {u -> RandomReal[]};  ivpsol = NDSolve[{     (* ODE system *)     {x''[t] == 0, y''[t] == -gravity},     (* initial values *)     {x[0] == y[0] == 0, x'[0] == vx0, y'[0] == vy0}     }, {x[t], y[t]}, {t, 0, duration}];  ParametricPlot[{x[t], y[t]} /. ivpsol, {t, 0, duration},   GridLines -> Automatic, AspectRatio -> 1/2]

Question: Now what options or language construct should i use to see approximations which are produced NDSolve while solving boundary-value-problem?

## Why can’t I solve a system of four third-order partial differential equations using NDSolve?

I solve the problem of optimal control. First, I solve a system of two equations (with respect to the variables p, q, V, n) with already optimal conditions (Kopt and Sopt) to check whether it is possible to solve such a system at all. It turns out that everything is OK, but the non-linear part of the equation cannot be discarded, it is important. I am using the NDSolve function with the StiffnessSwitching specification:

eq1 = D[n[x, t], t] + D[n[x, t]*V[x, t], x] == 0; eq2 = D[V[x, t], t] + Kopt[t]*x -      2*(2*250^2)*Sopt[t]*Sin[k*x]*Cos[k*x]*k + V[x, t]*D[V[x, t], x] +      g*D[n[x, t], x] -      D[ D[Sqrt[n[x, t]], x, x]/(2*Sqrt[n[x, t]]), x] == 0;  eqo1 = D[n[x, t], t] + D[n[x, t]*V[x, t], x] == 0; eqo2nl = D[V[x, t], t] + Kopt[t]*x -      8*(2*250^2)*Sopt[t]*Sin[k*x]*Cos[k*x]*k + V[x, t]*D[V[x, t], x] +      g*D[n[x, t],        x] - ((D[n[x, t], x])^3 -         2*n[x, t]*D[n[x, t], x]*D[n[x, t], {x, 2}] +         D[n[x, t], {x, 3}]*(n[x, t])^2)/(4*(n[x, t])^3) == 0;     cond1 = n[x, t] == ((Sqrt[Kzero]/Pi)^(1/2))*      Exp[-0.5*Sqrt[Kzero]*x^2] /. t -> 0; cond2 =  n[x, t] == 0. /. x -> -10; cond3 =  n[x, t] == 0. /. x -> 10;  cond4 = V[x, t] == 0. /. t -> 0; cond5 = V[x, t] == 0. /. x -> -10; cond6 = V[x, t] == 0. /. x -> 10;   sol5 = NDSolve[{eqo1, eqo2nl, cond1, cond2, cond3, cond4, cond5,      cond6}, {V, n}, {x, -10, 10}, {t, 0, Tv},     Method -> {"StiffnessSwitching",       Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 1,     PrecisionGoal -> 1];

Next, I solve four equations in order to find the control parameter K (depends on q), while Sopt is still known. The initial conditions are the same and as banal as possible. Add two more equations and K[t] formula:

K[t_] := -1*Integrate[q[x, t]*x^2, {x, -10, 10}];  eqo3 = D[q[x, t], t] == - n[x, t]*D[p[x, t], x]  -      V[x, t]*D[q[x, t], x];  eqo4 = D[p[x, t], t] +  D[p[x, t], x]*V[x, t] + g*D[q[x, t], x] == 0;

The linearized system has a solution, but not a nonlinear one! There are many problems, including the order of the equation:

The spatial derivative order of the PDE may not exceed two.

I try to remove the third order derivative from the nonlinear part:

The maximum derivative order of the nonlinear PDE coefficients for \ the Finite Element Method is larger than 1. It may help to rewrite \ the PDE in inactive form.

Why does he even use only FEM now? Why was a nonlinear system of two equations solved normally? What should I do, help me, please.

## Using NDSolve and PieceWise for boundary conditions for coupled PDEs

I’ve decided to try and use PieceWise in my boundary conditions for my 4 coupled PDEs, but it doesn’t seem to work out nicely. I’m receiving the error:

NDSolve::overdet: There are fewer dependent variables, {IPB[z,t],IPF[z,t],IS1B[z,t],IS1F[z,t]}, than equations, so the system is overdetermined.

Paramters:

ROCP = 1; RICP = 0; ROCS1 = 0.5; RICS1 = 1; L = 5 10^-2; n = 1.556; c = 3 10^10; lR = 3; R1 = 0.99;  gP = 20; gS1 = 20; sigmaSP = 10^-13;  toffset = 5 10^-9; twidth = 10 10^-9; Energy = 1.3 10^-13; rad = 0.015; PumpInt = ((Energy)/(twidth ))/(Pi rad^2) ; PumpPeak = Exp[0.5 ((-toffset)/(twidth))^2] PumpInt;  lc = n lR; alphaP = L/(2 lc); alphaS1 = (L - Log[R1])/(2 lc);  tmax = 20 10^-9;

The 4 PDEs I’m solving are:

PDE = {D[IPF[z, t], z] == -(n/c) D[IPF[z, t], t] -       gP IPF[z, t]*(IS1F[z, t] + IS1B[z, t]) - alphaP IPF[z, t],    -D[IPB[z, t], z] == -(n/c) D[IPB[z, t], t] -       gP IPB[z, t] (IS1F[z, t] + IS1B[z, t]) - alphaP IPB[z, t],    D[IS1F[z, t], z] == -(n/c) D[IS1F[z, t], t] +       gS1 IS1F[z, t] (IPF[z, t] + IPB[z, t]) - alphaS1 IS1F[z, t] +       sigmaSP (IPF[z, t] + IPB[z, t]),    -D[IS1B[z, t], z] == -(n/c) D[IS1B[z, t], t] +       gS1 IS1B[z, t] (IPF[z, t] + IPB[z, t]) - alphaS1 IS1B[z, t] +       sigmaSP (IPF[z, t] + IPB[z, t])};

And the Boundary conditions are:

BC = {IPF[z, t] ==      PieceWise[{{PumpPeak*Exp[-0.5 ((t - toffset)/twidth)^2],         z == 0 && t > 0}, {0, 0 <= z <= lR && t == 0}}],    IPB[z, 0] == 0, IPB[lR, t] == IPF[lR, t] ROCP,    IS1F[0, t] == IS1B[0, t] RICS1, IS1B[lR, t] == IS1F[lR, t] ROCS1,     IS1F[z, 0] ==  0, IS1B[z, 0] ==  0};

And the NDSolve:

solInt = NDSolve[{PDE, BC}, {IPF, IPB, IS1F, IS1B}, {z, 0, lR}, {t, 0,      tmax}, MaxStepFraction -> {1/750, 1/10},     Method -> {"BDF", "MaxDifferenceOrder" -> 5},     InterpolationOrder -> All];

The main idea behind the boundary conditions is that IPF[z,t] = 0 at time t = 0 and for any z in the range 0 – lR, but for t > 0 and z = 0, it will have the form of the Gaussian equation.

Is it possible to impose these initial conditions?

Side note: Sorry if my post are similar, I didn’t know if I should piggyback on my last question or if this was something I should ask from scratch.

Cheers!

## How to adjust PrecisionGoal and AccuracyGoal to solve my equation using NDSolve[]?

This is a very interesting question. I’m solving a differential equation to get a function z[$$\rho$$], the expected solution needs to start from $$z[\rho_c]=z[\frac{13}{5}]=\frac{24195890215702774518563237183870329}{8112963841460668169578900514406400}\approx 2.98$$, until when the function $$z[\rho]$$ hits on the $$\rho$$-axis. The code is as follows:

precTmp = 4 MachinePrecision; sol = (\[Rho]c = 13/5;    vt = Sqrt[zt^2 - \[Rho]c^2] - zt;     ParametricNDSolveValue[{(1 - z[\[Rho]]^3)^2/z[\[Rho]] +         3/2 z[\[Rho]]^2 Derivative[1][          z][\[Rho]]^2 + (1 - z[\[Rho]]^3) (Derivative[1][z][\[Rho]]^2/           z[\[Rho]] + (z^\[Prime]\[Prime])[\[Rho]]) + (        169 (-(13/5) + zt)^2 (13/5 +            zt)^2 z[\[Rho]]^2 (-3 z[\[Rho]]^2 +            2 (z^\[Prime]\[Prime])[\[Rho]]))/(200 zt^2) == 0,       Derivative[1][z][13/5] == -((13 (2 - (-(169/25) + zt^2)^(3/2)))/(        10 Sqrt[-(169/25) + zt^2])), z[13/5] == Sqrt[-(169/25) + zt^2],       WhenEvent[z[\[Rho]] == 10^-3, "StopIntegration"],       WhenEvent[z[\[Rho]] == 1.1 Sqrt[zt^2 - \[Rho]c^2],        "StopIntegration"]}, z, {\[Rho], \[Rho]c, R + 1}, {zt},      WorkingPrecision -> precTmp, AccuracyGoal -> precTmp/2, MaxSteps -> 10^6,      Method -> "StiffnessSwitching"]);

Note that we can adjust the parameter precTmp to ramp up the precision. To check the correctness of the solution, it’s natural to plot the solution out, as follows:

ztValue = 24195890215702774518563237183870329/   8112963841460668169578900514406400; Plot[sol[ztValue][\[Rho]], {\[Rho], sol[ztValue]["Domain"][[1, 1]],    sol[ztValue]["Domain"][[1, 2]]}, PlotRange -> All,   WorkingPrecision -> precTmp]

The first thing I can’t understand is that the shape of the solution $$z[\rho]$$ is different when I change the precision, for example, taking precTmp to be 3 MachinePrecision, 5 MachinePrecision, or 6 MachinePrecision, the solution is decreasing and hit on $$\rho$$-axis, but when taking precTmp to be 4 MachinePrecision, 8 MachinePrecision or 12 MachinePrecision, the solution turns out to be increasing. My first question is which one of them is the correct answer?

The second question that I can’t understand is when providing a PrecisionGoal in ParametricNDSolveValue[] function, for example, we take precTmp = 4 MachinePrecision again, and add PrecisionGoal->precTmp, the solution becomes a decreasing function, which is originally an increasing one. I’ve already learned the difference between precision and accuracy, which can somehow be regarded to be equivalent as long as the function doesn’t go to zero. Thus my second question is how to understand the difference brought by the PrecisionGoal set here.

## Error while giving initial condition in NDSolve

I am trying to solve following differential equation $$\frac{d^2y}{dx^2}+(a+b \frac{2}{\pi}\tan^{-1} x)y=0$$ with the initial condition $$y(-10)=e^{10i\sqrt{a+b}}$$ and $$y'(-10)=-i\sqrt{a+b}e^{10i\sqrt{a+b}}$$.

To implement it I write the following code

s = ParametricNDSolve[{y''[x] + (a + b (2 + 2/Pi ArcTan[x])) y[x] ==  0, y[-10] = Exp[I 10 Sqrt[a + b]], y'[-10] = -I Sqrt[ a + b]*Exp[ I 10 Sqrt[a + b]] }, y, {x, -10, 10}, {a, b}]

But I am getting an error that ParamatericNDSolve expects equation or list of equations instead of $$e^{10i\sqrt{a+b}}$$ in the first argument. Can anyone point me out where am I making the mistake?

## Trouble implementing outflow boundary condition when trying to solve a pde using NDSolve

Trying to solve the following pde: $$\partial_{t}y + c\partial_{c}y = 0$$ (for simplicity $$c=1$$).

For the initial data I am using a Gaussian. The problem surges when I am trying to implement the outflow boundary condition as it was suggested to me, namely $$\frac{\partial{y}}{\partial{t}} =0$$ at $$x =0$$.

So far my code is pretty simple:

v = 1 ; L = 2;   With[{y = y[t, x]},   eq = D[y, t] + vD[y, x] == 0;  ic = y == Exp[-x^2] /. t -> 0;  bc = {D[y, x] == 0 /. x -> 0 }];    mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines",   "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,   "MinPoints" -> n, "DifferenceOrder" -> o}};   sol = NDSolveValue[{eq, ic, bc}, y, {t, 0, 1}, {x, 0, L},   Method -> mol[100, 4]];  {t0, tend} = sol["Domain"][[1]];   Manipulate[  Plot[sol[t, x], {x, 0, L}, PlotRange -> {-10, 10}], {t, 0, tend}];

It stems from answers to questions previously asked here, and I intend later to test finite difference methods (BTCS/FTCS) as done in Schemes for nonlinear advection equation.

However I am not being able to evolve the equation do to confusion when trying to implement the BC, I get the following:

NDSolveValue: Boundary condition $$y^{(0,1)}[t,0]$$ should have derivatives of order lower than the differential order of the partial    differential equation.

This is expected as I am not sure what would be the best way to impose BC on the problem.

If anyone has any suggestions they would be welcolmed.

Thanks.

## Save NDSolve output before break down point?

When NDSolve breaks down at a certain point and returns

NDSolve::ndstf: At t == 0.000292874155196631662368876027253458083949302914788253674940., system appears to be stiff. Methods Automatic, BDF, or StiffnessSwitching may be more appropriate.

is there a way to save or plot the solution from t=0 up to t == 0.0002 and see what went wrong?

## How to NDSolve this equation?

I am trying to find the answer to f[r] in the following differential equation. Is there anything wrong? Why the output is printed with error?

eq=-2 + r - r*f[r] +    (0.01*(-2*Derivative[1][f][r]*(6*f[r] - 6*f[r]^2 +         r*Derivative[1][f][r]*(3 + r*Derivative[1][f][            r])) + 6*r*f[r]*(2 - 2*f[r] +         r*Derivative[1][f][r])*Derivative[2][f][r]))/    r^2;  NDSolve[{eq==0},f,{r,0,5}]

All the additional informations are $$f(0)=\infty$$, and $$f(\infty)=1$$.

## Export NDSolve output to .txt

I tried to store the output of the following piece of code

s = NDSolve[{u'[x] ==  v[x] , v'[x] == u[x], u[0] == 1,      v[0] == 1}, {u, v}, {x, 0, 1}, WorkingPrecision -> 70,     PrecisionGoal -> 20] ;

in a .txt file with the commands

SetDirectory[NotebookDirectory[]]; data = Table[{x, u[x] /. s, v[x] /. s}, {x, 0, 1, 0.1}] // TableForm Export["data.txt", data, "List"];

Unfortunately the resulting .txt had the following ugly arrangement

0.{1.}{1.} 0.1{1.1051709180440972}{1.1051709180440972} 0.2{1.2214027580171631}{1.2214027580171631} 0.30000000000000004{1.349858807431976}{1.349858807431976} 0.4{1.4918246976372467}{1.4918246976372467} 0.5{1.648721270402639}{1.648721270402639} 0.6000000000000001{1.8221187976152273}{1.8221187976152273} 0.7000000000000001{2.013752707127757}{2.013752707127757} 0.8{2.2255409261449164}{2.2255409261449164} 0.9{2.4596031083207683}{2.4596031083207683} 1.{2.718281828459045}{2.718281828459045}

I want to get something like this

0.0    1.00000  1.00000 0.1    1.10517  1.10517 0.2    1.22140  1.22140  0.3    1.34985  1.34985  ...

Can anyone help?