Second order boundary value problem numerical solution

Hello,

I have a second order boundary value problem with drichlet boundary conditions as follows:

(d^2(m(y))/dy^2) = m(y)*(c1+c1*P/b(y))

domain [0 H];

Boundary Conditions

m(0) = 0;

m(H) = 0;

In this case b(y) is a known vector; it has a (1Xy_mesh) dimensions ; The domain of b(y) and m(y) are the same. Only thing is that b(y) is an input for solving m(y); c1 and P are constants

How can i solve this equation?

Answer 1

[Riccardo Scorretti]

Hi. There are plenty of methods to solve your problem. Perhaps, the simplest one is by using Finite Difference; you can find many excellent textbooks on this topics. In particular, I suggest Numerical approximation of partial differential equations by by Quarteroni and Valli.

That's being said, with these boundary conditions the analytical solution is because: , and obviously hence: .

Assuming that indeed you want to solve: here is a small code which executes the computation. The idea is to approximate: , where , = nodes of the compuational domain, and = discretization step.

Hence you can discretize your PDE, and this boils up into a linear system:

together with imposed Dirichlet boundary conditions:

% Parameters of the PDE
% ---------------------
c1 = 1;
H = 2;
P = 1;
b = @(y)  1+y;        % Source term (it has to be different from 0)
m0 = 0;               % Value of m(0)
mH = 0;               % Value of m(H)

% Parameters of the solver
% ------------------------
N = 100;                % Number of nodes

% Finite Difference solver
% ------------------------
A = sparse(N, N);       % Matrix of the linear system
z = zeros(N, 1);        % Known vector 

h = H/(N-1);            % = Discretization step
y = linspace(0, H, N);  % = Coordinates of the nodes
bn = feval(b, y);       % = Source term b(y) computed in each node

% Assembly the laplacian operator
one = ones(N, 1);
A = spdiags([one -2*one one]./h^2, -1:1, N, N);

% Assembly the term -c1*m
A = A + spdiags(-c1*one, 0, N, N);

% Assembly the source term
z = c1*P./bn(:);

% Impose Dirichlet boundary conditions by replacing the first and last
% equations of the linear system
A(1,:) = 0 ; A(1,1) = 1 ; z(1) = m0;     % impose m(1) = 0
A(N,:) = 0 ; A(N,N) = 1 ; z(N) = mH;     % impose m(N) = 0

% Solve the linear system and plot the result
m = A \ z;
figure
plot(y, m, '.-');
xlabel('y') ; ylabel('m(y)');

You can modify the code so as to check it with more trivial problems, for instance: with and , which will provide the (expected) solution .

Enjoy.

[Oguz Altunkas]

Thank you for your answer,

How did you decide b(y) to be equal to 1+y. Because, we have a vector of numbers at each node. We do not know the b(y) vector in terms of y. I am confused at this point?

Thank you in advance

[Riccardo Scorretti]

Well, it was just an example. If you already have the values of b in each node, of course you can just use it place of the variable bn (and use your own mesh instead of variable y).

If you know the value of b for a set of points, and you want to solve the problem with a finer discretization, you will have to interpolate the value of b on the nodes (= variable y).

This kind of details depend specifically on your program.