Boundary conditions with stiff problems

Hello, I'm trying to solve a system of PDEs that depends on space (1 dimension) and time. I want to solve it using a second order discretization in space and then using ode23s to solve the system of ODEs (method of lines).

The problem is that I have a laplacian operator and I can just integrate in the interior nodes (let's say from 2 to n-1). The values for nodes 1 and n depend on their neighbours and should be updated at each time step.

How could I set these restrictions to solver?

Here is the scheme:

function dydt = fun(t,u)

for j = 2:n-1
impose right hand side function
end

%Now I want to impose the value in y(1) and y(end)

end

Answer 1

[Torsten]

Boundary conditions don't depend on neighbour values, but are given independently.

What are the boundary conditions for your PDE ?

Best wishes

Torsten.

[Albert Jimenez]

The boundary conditions are du/dx = G(u), for a given function G at x=0 and x=1 (the same for two points).

I have tried to discretize them too but the solution is not as expected (in fact, I obtain a warning: matrix is singular, close to singular or badly scaled).

I will send you the paper I'm trying to simulate if needed.

Thanks

Answer 2

[Torsten]

Use the equations

(u(2)-u(1))/(x(2)-x(1)) = G(u(1))
(u(n)-u(n-1))/(x(n)-x(n-1)) = G(u(n))

to solve for u(1) (u at 0) and u(n) (u at 1).

Then use these values for the discretization in the inner grid points.

Best wishes

Torsten.

Answer 3

[Albert Jimenez]

The file attached contains the system I want to solve (1D). Using a second order discretization in space, the laplacian looks like

(u(j-1)+u(j+1)-2u(j))/(deltax^2)

The boundary conditions are

du/dx = A(u-B) at x=0 and 1

where A and B are constants.

If I discretize the first order derivative with a second order discretization, I obtain

(u(0)-u(2))/(2*deltax) = A(u(1)-B)

where u(0) is the value of u in a ficticious node. Then I can isolate the value of u(0) and replace in the scheme equation for the first node.

Is this what you were suggesting? This way, and assuming I don't have any mistake in the code, I get constant solutions (and I shouldn't).

Here is the code I am using: http://pastebin.com/ANe9ZgUz

Thanks

[Torsten]

If your boundary condition reads

du/dx = A*(u-B) at x=0 and x=1,

then the discretized form is

(u(2)-u(0))/(2*deltax) = A*(u(1)-B)
(u(n+1)-u(n-1))/(2*deltax) = A*(u(n)-B)

So you have a sign error in your equation from above.

Best wishes

Torsten.

[Albert Jimenez]

Yes! You are right. In fact, I could express the value on the first and last node in terms of the second and third this way:

u(1) = (4u(2)-u(3)+C1)/C2

where C1 and C2 depend on A and B. Analogously for u(n). However, this does not solve my problem.

I have tried to integrate in an small interval and then impose these values at the end, but Matlab says that the matrix is singular, close to singular or badly scaled (don't know what matrix, it seems to be something for solving a linear system inside ode23s function).

Thanks. I will comment if I find the solution.

[Torsten]

Your loop

for j=2:n-1

must run from 1 to n.

Best wishes

Torsten.

[Albert Jimenez]

Don't think so... I have a Laplacian and I would point to values u(-1) and u(n+1). The boundary points should be integrated out of the loop.

[Torsten]

No. The boundary conditions give you u(0) depending on u(1),u(2) and u(n+1) depending on u(n-1),u(n).

Inserting these expressions in the discretized Laplacian in x(1) and x(n) give you the time derivatives for u(1) and u(n).

Best wishes

Torsten.

[Albert Jimenez]

Thanks! You are right and this is what I tried at the begining. However, I do not obtain the result as shown in the paper I am trying to reproduce. I could sent you the paper if you are interested.

Thanks again.