4
$\begingroup$

I have two coupled differential equations as follows $$ \frac{\partial }{\partial x}U(x,y) =2V(x,y), $$ $$ \frac{\partial }{\partial y}V(x,y) =V(x,y)U(x,y)+1, $$

with conditions $$ U(0,y)=sin(y),V(x,0)=0. $$

We want to solve these equations numerically with mathematics. I used the NDSolve command as follows

NDSolve[{D[U[x, y], x] == 2 V[x, y],D[V[x, y], y] == V[x, y] U[x, y] + 1, V[x, 0] == 0,U[0, y] == Sin[y]}, {U, V}, {x, 0, 2}, {y, 0, 2}]

Problem

When I run this, I recieve the following errors:

NDSolve`FEM`InitializePDECoefficients::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.
FindRoot::stfail: The method AffineCovariantNewton failed to compute the next step.
FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option.

Can you help me fix this?

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ If you look at the first error message on your screen, you'll see a blue 'i' icion at the right. Click on that and it takes you to the error message description: FEMDocumentation/ref/message/InitializePDECoefficients/femcscd. If you scroll t the bottom, there is a recommendation to use a finer mesh, and an example of how to do so. The other FindRoot errors come later, probably as a consequence of the first two. So maybe try setting up a fine mesh x-y region, and see if that helps at all. $\endgroup$ Commented Jul 17, 2023 at 2:57
  • $\begingroup$ Also, at the very top of the page, there's a link to "FEM Best Practices" which describes how to set up artificial diffusion, as suggested by the message. the example shows a set of first order equations for rho[x,t] and v[x,t], similar to yours (though without the nonlinear term). What is done is to add a small second order spatial derivative term (const rho[t,x],{x,2}] etc.) to each equation, adjusting the constants to smooth out the solution but keeping the extra terms as small as possible. $\endgroup$ Commented Jul 17, 2023 at 3:35
  • $\begingroup$ @A.D Where do the equations come from? $\endgroup$ Commented Jul 19, 2023 at 8:19

1 Answer 1

Reset to default
1
$\begingroup$

I tried adding a pair of "artificial diffusion" terms, as described in my comment above. Your case differs from the example in FEM Best Practices. The exmaple uses two time derivatives, you use an x derivative for u and a y derivative for v. So, I'm not sure what I did is correct, but you can modify the code and try out different combinations to see what works and what doesn't.

I changed your variables to lower case, since a number of upper case single letter variables are reserved symbols in Mathematica. Not U,V as far as I know, but it's a good habit to avoid single letter upper case variable names in general.

The first two equations are artificial diffusion terms for u,v. 'cdif' is a constant to be varied. Add each term to the corresponding pde term, pick a value for cdif, and try running NDSolve. My initial choice was cdif = .01, but that failed inside FindRoot due to convergence difficulties. You might try setting up a plane region with a fine mesh to see if that helps.

adu = cdif D[u[x,y],{y,2}];

adv = cdif D[v[x,y],{x,2}];

pde={D[u[x,y],x]-2 v[x,y]==adu,D[v[x,y],y]-(v[x,y] u[x,y]+1)==adv};

bc = {v[x,0]==0,u[0,y]==Sin[y]};

This didn't generate warnings. I'm not going to copy the InterpolatingFunctions over - they end up in OutputForm and are huge.

cdif=.1;NDSolve[Flatten[{pde,bc}],{u,v},{x,0,2},{y,0,2}]

Here are plots of u and v:

3d plot of u[x,y] 3d plot of v[x,y]

I have no idea if this is what your were expecting, but you have something to tweak.

EDIT: I tried experimenting with a finer mesh, and using different artificial diffusion terms, such as Div[Grad[u]], etc., but saw no improvement; if anything, it made things worse. I did find that cdif can be lowered to roughly 0.075 before DSolve fails, but the plots aren't too different.

The artificial diffusion approach is a somewhat standard method used in computational fluid dynamics. If your equations model a fluid, then you can justify adding such terms on physical grounds since a realistic fluid would exhibit such behavior. Since you're modifying the equations, you want to keep the added terms as small as possible to avoid corrupting your equations too much.

With that in mind, I took the solutions with cdif = 0.075 and made some plots to compare the relative sizes of various terms in the equation. It's not a comprehensive analysis by any means. I arbitrarily took a value of x -> 0.5, and plotted against y over the whole range.(1d plots are easier to read, at least for me.)

Plot[Evaluate[{u[x,y],v[x,y]}/.Flatten@nds/.x->.5],{y,0,2},PlotLegends->{"u","v"}]

u,v plot

Plot[Evaluate[{D[u[x,y],x],D[v[x,y],y]}/.Flatten@nds/.x->.5],{y,0,2},PlotLegends->{"du/dx","dv/dy"}]

u',v' plot

Plot[Evaluate[{cdif D[u[x,y],{y,2}],cdif D[v[x,y],{x,2}]}/.Flatten@nds/.x->.5],{y,0,2},PlotLegends->{"d^2u/dy^2","d@v/dx^2"}]

u'',v'' plot

The very tentative conclusion I draw is that the extra diffusion terms aren't small compared to u,v, but are small compared to the first derivatives. Unfortunately, as mentioned, I had no luck in shrinking them down further, without having NDSolve exit with no solution.

Improve this answer
$\endgroup$
5
  • $\begingroup$ It seems unusual to me to change a first order pde-system to a second order system with small parameter cdif. It looks like a singular system? $\endgroup$ Commented Jul 17, 2023 at 10:12
  • 1
    $\begingroup$ @Ulrich Neumann It's numerically unstable, but not necessarily singular. You get a solution, but there's a big ridge on one boundary for u. I'm not very knowledgeable about this, but following the link from the error msg to the documentation, this (and using a finer mesh) seems to be the recommended solution.The error message suggests adding artificial diffusion, and adding a second order term is what's shown as the suggested approach. When I have time, I'll try to research this further and try to understand it better. $\endgroup$ Commented Jul 17, 2023 at 15:30
  • $\begingroup$ Thanks for your reply, I'm very curious about the result $\endgroup$ Commented Jul 18, 2023 at 7:25
  • $\begingroup$ Interesting edit! I tried to solve with a refined boundary mesh, but solution looks totally different $\endgroup$ Commented Jul 19, 2023 at 8:22
  • $\begingroup$ @UlrichNeumann I just tried a simple rectangular mesh with finer spacing, but NDSolve failed to get anything. I've run out of ideas for the time being. If you're interested in this problem, maybe you'll find a way to improve on this. Good luck! $\endgroup$ Commented Jul 19, 2023 at 16:32

Your Answer

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.