Solving the Implicit Euler ODE with Boundary Conditions

[member 428835]

Homework Statement

Write an implicit Euler code to solve the system $c'(x) = \epsilon c''(x)-kc(x)$ subject to $1-c(0)+\epsilon c'(0) = 0$ and $c'(1)=0$.

Homework Equations

Nothing out of the ordinary comes to mind.

The Attempt at a Solution

In the following code, there is central difference implicit Euler being used, taking the form $A*Y=B$ where $Y = [c_0,c_1...c_{N+1}]^T$. I understand the procedure except when considering the boundary conditions. I marked the three lines corresponding to this confusion with a HOW TO ARRIVE comment. This isn't the entire script, but this is the part I need help with. There is the BC at $x=1$ but I think if I understand how the $x=0$ BC works I can figure out $x=1$.

A=zeros(nstep+1,nstep+1);           % initialize A matrix (A*Y=B)
B=zeros(nstep+1,1);                 % initialize B matrix

% components of the matrix A (A*Y=B)
a1=eps/h^2+1/(2*h);
b1=-2*eps/h^2-k;
c1=eps/h^2-1/(2*h);

% configuration of the matrix A
A(1,1)=b1-a1*2*h/eps;% HOW TO ARRIVE
A(1,2)=c1+a1;% HOW TO ARRIVE
B(1,1)=-2*a1*h/eps;% HOW TO ARRIVE
for i=2:N
    A(i,i-1)=a1;
    A(i,i)=b1;
    A(i,i+1)=c1;
end
A(N+1,N)=a1+c1;
A(N+1,N+1)=b1;

AI=inv(A);                      % matrix inversion by Matlab func.
Y=AI*B;                         % final solution of y

 

[DrClaude]

To remove some confusion, let's take instead $Y = Y = [y_1,y_2...y_{N+1}]^T$ to have the same index notation as the matrix and to keep $c_1$ for the variable $\texttt{c1}$.

Away from the boundary, you have that element $i$ of $AY$ is $a_1 y_{i-1} + b_1 y_i + c_1 y_{i+1}$. Find out what you get for the first element of $AY$ and compare.

 

[Chestermiller]

When you say you are using implicit Euler, does that mean you are solving this problem iteratively? Let's see your difference equations.

 

[member 428835]

Away from the boundary, you have that element $i$ of $AY$ is $a_1 y_{i-1} + b_1 y_i + c_1 y_{i+1}$. Find out what you get for the first element of $AY$ and compare.

The first element would satisfy $a_1 y_{1} + b_1 y_2 + c_1 y_{3} = 0$ and also the boundary condition $1-y(0)+\epsilon y'(0)$, which is expressed $1 - y_1+\epsilon (y_3-y_1)/dx = 0$. Substituting this BC into the difference equation yields $y_1(a_1+c_1+c_1 dx/\epsilon)+b_1y_2 = c_1 dx/\epsilon$. Then shouldn't $A_{11} = a_1+c_1+c_1 dx/\epsilon, \: A_{12} = b_1, \: A_{13} = 0$ and $B_{11} = c_1 dx/\epsilon$?

When you say you are using implicit Euler, does that mean you are solving this problem iteratively? Let's see your difference equations.

It's here

% dy/dx caculation with central difference scheme
Y0=2*h/eps-2*h/eps*Y(1,1)+Y(2,1);
DY(1,1)=(Y(2,1)-Y0)/(2*h);
for i=2:N                       
    DY(i,1)=(Y(i+1,1)-Y(i-1,1))/(2*h);
end
YN2=Y(N,1);
DY(N+1,1)=(YN2-Y(N,1))/(2*h)

 

[Chestermiller]

The first element would satisfy $a_1 y_{1} + b_1 y_2 + c_1 y_{3} = 0$ and also the boundary condition $1-y(0)+\epsilon y'(0)$, which is expressed $1 - y_1+\epsilon (y_3-y_1)/dx = 0$. Substituting this BC into the difference equation yields $y_1(a_1+c_1+c_1 dx/\epsilon)+b_1y_2 = c_1 dx/\epsilon$. Then shouldn't $A_{11} = a_1+c_1+c_1 dx/\epsilon, \: A_{12} = b_1, \: A_{13} = 0$ and $B_{11} = c_1 dx/\epsilon$?
It's here

% dy/dx caculation with central difference scheme
Y0=2*h/eps-2*h/eps*Y(1,1)+Y(2,1);
DY(1,1)=(Y(2,1)-Y0)/(2*h);
for i=2:N                      
    DY(i,1)=(Y(i+1,1)-Y(i-1,1))/(2*h);
end
YN2=Y(N,1);
DY(N+1,1)=(YN2-Y(N,1))/(2*h)

So, you are solving this split boundary value problem using the "shooting method?"

 

[member 428835]

So, you are solving this split boundary value problem using the "shooting method?"

What's a split bvp? And yep!

 

[Chestermiller]

You seem to be using a 2nd order central difference finite difference approximation when Euler is a first order difference.

 

[member 428835]

You seem to be using a 2nd order central difference finite difference approximation when Euler is a first order difference.

So you're saying I'm inconsistent on some level?

 

[Chestermiller]

So you're saying I'm inconsistent on some level?

I guess yes. Euler is a one-sided difference approximation. Do you want to solve this problem with central difference or with Euler?

 

[member 428835]

I guess yes. Euler is a one-sided difference approximation. Do you want to solve this problem with central difference or with Euler?

Sorry it took me so long to reply! How about both?

 

[Chestermiller]

Let's make sure we're on the same page. You're solving a coupled set of linear ODEs of the form:$$\frac{dy}{dx}=Ay$$ subject to boundary conditions at x=0 and x=1, where y is the solution vector and A is a real square matrix. You are using the shooting method to satisfy the BC at y = 1, and using implicit Euler for the differential equation: $$y^{n+1}-\Delta xAy^{n+1}=y^n$$

Are we in agreement?

 

[member 428835]

Let's make sure we're on the same page. You're solving a coupled set of linear ODEs of the form:$$\frac{dy}{dx}=Ay$$ subject to boundary conditions at x=0 and x=1, where y is the solution vector and A is a real square matrix. You are using the shooting method to satisfy the BC at y = 1, and using implicit Euler for the differential equation: $$y^{n+1}-\Delta xAy^{n+1}=y^n$$

Are we in agreement?

Yes we are :)

 

[Chestermiller]

OK. So we set $y_1=c$ and $y_2=\frac{dc}{dx}$. The our coupled differential equations become: $$\frac{dy_2}{dx}=\frac{1}{\epsilon}y_2+\frac{k}{\epsilon}y_1$$and$$\frac{dy_1}{dx}=y_2$$This is subject to the boundary conditions:
$$y_1(0)=1+\epsilon y_2(0)$$and$$y_2(1)=0$$. Applying the implicit Euler relationship to the differential equations, we obtain:

$$y_1^{n+1}-\Delta xy_2^{n+1}=y_1^n$$
$$-\Delta x\frac{k}{\epsilon}y_1^{n+1}+(1-\Delta x\frac{1}{\epsilon})y_2^{n+1}=y_2^n$$
These represent two linear algebraic equations in the two unknowns for $y_1^{n+1}$ and $y_2^{n+1}$ at grid point n+1 in terms of $y_1^{n}$ and $y_2^{n}$ to grid point n.

OK so far?

 

[member 428835]

This is subject to the boundary conditions:
$$y_1(0)=1-\epsilon y_2(0)$$and$$y_2(1)=0$$

I think the first BC is $$y_1(0)=1+\epsilon y_2(0).$$ Otherwise I am good!

 

[Chestermiller]

I think the first BC is $$y_1(0)=1+\epsilon y_2(0).$$ Otherwise I am good!

Oops. Sorry. I'll go back and change that.

Okay, so you can use these two equations over and over again from one boundary to the other to determine the values of the solution vector components at each successive grid point from those at each previous grid point.

OK with this?

We can discuss the implementation of the shooting method once we have agreement up to this point.

 

[member 428835]

Oops. Sorry. I'll go back and change that.

Okay, so you can use these two equations over and over again from one boundary to the other to determine the values of the solution vector components at each successive grid point from those at each previous grid point.

OK with this?

We can discuss the implementation of the shooting method once we have agreement up to this point.

Yea, I agree with this. But the second BC hasn't been addressed. The shooting method would imply we guess a value for $y_1(0)$, which then determines $y_2(0)$ and see how this agrees with $y_2(1) = 0$, right?

 

[Chestermiller]

Yea, I agree with this. But the second BC hasn't been addressed. The shooting method would imply we guess a value for $y_1(0)$, which then determines $y_2(0)$ and see how this agrees with $y_2(1) = 0$, right?

Right, except, for some reason, I would think of it as guessing the value for $y_2(0)$, which determines $y_1(0)$. The initial guess I would use for $y_2(0)$ would be zero. Would you know what to do next, once you determined $y_2(1)$, for the guessed value at x = 0?

 

[member 428835]

Would you know what to do next, once you determined $y_2(1)$, for the guessed value at x = 0?

Since $y_1 = c$ and from BC we have determined $y_1$ is there any more work to do?

 

[Chestermiller]

Since $y_1 = c$ and from BC we have determined $y_1$ is there any more work to do?

I think so. It starts off with $y_1^{(0)}=1$, and $y_2^{(0)}=0$. But then you get $y_1^{(1)}=1$ and $y_2^{(1)}\neq 0$ at $x = \Delta x$ . Then the integration continues up to x = 1.

 

[member 428835]

I think so. It starts off with $y_1^{(0)}=1$, and $y_2^{(0)}=0$. But then you get $y_1^{(1)}=1$ and $y_2^{(1)}\neq 0$ at $x = \Delta x$ . Then the integration continues up to x = 1.

Perhaps I'm not understanding you, but am I right in thinking it doesn't matter what value $y_2^{(1)}$ is so long as $y_2^N = 0 : y_2^N \equiv y_2(1)$. If this does not happen we guess a new value for $y_1^{(0)}$ and reiterate, right?

 

[Chestermiller]

Perhaps I'm not understanding you, but am I right in thinking it doesn't matter what value $y_2^{(1)}$ is so long as $y_2^N = 0 : y_2^N \equiv y_2(1)$. If this does not happen we guess a new value for $y_1^{(0)}$ and reiterate, right?

I think we agree. It's just that I would be iterating on the guess for $y_2(0)$, not on a guess for $y_1(0)$. I realize that the two approaches are equivalent, but my intuition tells me to focus on $y_2(0)$, not $y_1(0)$. In doing the integration, I would also make sure that I chose $\Delta x$ such that $\Delta x\lt \epsilon$; otherwise, the numerical solution is going to have oscillations.

 

[member 428835]

In doing the integration, I would also make sure that I chose $\Delta x$ such that $\Delta x\lt \epsilon$; otherwise, the numerical solution is going to have oscillations.

Can you explain how you know this? Also, I assume you select $y_2$ as the guess since we are given $y_2$ at the far boundary?

 

[Chestermiller]

I know this because of the coefficient of y2 in the 2nd difference equation. I select y2 as the guess because it is multiplied by epsilon (which presumably is small).

The original differential equation has all the characteristics of a singular perturbations system, with epsilon multiplying the highest order derivative. I would expect the solution to be varying very rapidly near one of the boundaries if epsilon is small.

 

[member 428835]

Thanks a lot!

 

[Chestermiller]

Thanks a lot!

Try solving analytically for the case where epsilon is zero and k is very large to get you started (to see what to expect).