Is the oscillation in the picture a boundary condition problem?

Hi,

I have an equation which is as follows

When I solve this equation using the pdepe it shows an oscillation pattern at the edge as it is here (shown by arrow).

I changed my parameters a lot but this oscillation sometime reduces, but not zero. My bc are as follows

pl =  0;
ql =  1;

pr = 0;
qr = 1;

Can you tell how to reduce this oscillations?

Also, in the above equation, for the diffusion coefficient (sigma) the minimum value for which the pdepe converges is 1 below which the solution becomes diverging, it should not be actually.

with regards

[Aquatris]

Provide the full code so it is easier to analyze.

[Torsten]

Which boundary condition do you want to set at the right endpoint ? You set I * dI/dx = 0 - this can mean I = 0 or dI/dx = 0.

[Rahul]

Hi Torsten,

I want only dI/dx to be =0 at the right boundary but I should not be zero.

[Torsten]

Then MATLAB chose the wrong possible option ... I don't think you can decide whether you want to set I = 0 or dI/dx = 0 if you choose p = 0 and q = 1. But I also don't see a way to specify dI/dx = 0 uniquely with "pdepe".

[Rahul]

Hi Torsten,

Thanks a lot for your reply.

I'd like to know, in case I want to keep both I and dI/dx to be not equal to zero,in that case how should I write the bc for right boundary?

[John D'Errico]

I'll suggest you will never kill the oscillations off completely. However, you can probably reduce them I think, merely by use of a finer xmesh and/or tspan.

[Torsten]

I'd like to know, in case I want to keep both I and dI/dx to be not equal to zero,in that case how should I write the bc for right boundary?

You could choose a "mixed" condition - thus a condition where both p and q are not equal to 0.

But can you show us your code ? Maybe there is some other reason for the solution's behaviour.

[Rahul]

Hi Torsten,

I'd like to share the code to you in your email and not in public.

[Torsten]

Sorry, in this case it's not possible. I don't want to share my email.

Just to be sure that pr = 0 and qr = 1 imply what I meant it would imply:

You set

c = 1
f = sigma*I*dIdx
s = chi*(mu*I-gamma_NL*I^2)

in pdefun, didn't you ?

[Rahul]

Hi Torsten,

I'm using it exactly.

But I have a query. Since the equation is a convection-diffusion equation, if I wish to use Robin bc, how should I proceed, can you please advice me?

Thanks.

[Umar]

Hi @Rahul,

As I understand that you are probably thinking that reason for your oscillations are due to arising from insufficient resolution or inappropriate boundary conditions, that is why you keep asking @Torsten to suggest ideas to resolve this issue by asking him “ case I want to keep both I and dI/dx to be not equal to zero,in that case how should I write the bc for right boundary?” and round Robin, so for your boundary conditions, if you want to maintain both the variable ( I ) and its derivative dI/dx non-zero at the right endpoint, you can proceed implementing Robin boundary conditions as follows

% Left and right   boundary conditions

bc = @(xl, ul, xr, ur) deal(ul - 0, qr*ur - pr);

So, in this code snippet, qr and pr are coefficients that you can adjust to control the behavior of the solution at the boundary. I will also suggest that diffusion coefficient sigma is set to a minimum of 1, as values below this can lead to divergence. Additionally, consider using a higher-order method or applying artificial viscosity to dampen oscillations. This approach will help stabilize your solution while maintaining accuracy.

Afterwards, you mentioned “the equation is a convection-diffusion equation, if I wish to use Robin bc, how should I proceed,”, I would proceed by discretizing the domain by defining spatial grid and time steps. Then, setup the equations for interior points by using the standard finite difference approximations for convection and diffusion. Then, at the boundary, replacing the standard boundary conditions with the Robin condition. For example, at the left boundary (i=1):

u(1) = (g - beta * (u(2) - u(1))/dx) / alpha;

Finally, updating your solution iteratively, so that the Robin condition is applied at each time step.

[Torsten]

@Rahul

Differentiating your second-order term gives I*d^2I/x^2 + (dI/dx)^2. So you should avoid I to become zero - this would introduce a singularity into your problem.

Maybe you should test how the steady-state equation (dI/dt = 0) behaves. Try bvp4c.

[Rahul]

Hi,

Thankyou everyone. You all helped me a lot by providing me valuable insights for the problem I encountered.

rgds

[Umar]

Hi @Rahul,

Thank you for your kind words. I am glad to hear that the insights shared were helpful in addressing the problem you encountered. It is always a pleasure to collaborate and support one another. If you have any further questions or need additional assistance, please do not hesitate to reach out.

[Umar]

Hi @Rahul,

The code that you provided was incorrect. So, I had to fix the errors and update the code. Please see attached. Hopefully, this answers your question.

[Rahul]

Hi Umar,

Thanks a lot for your support.