Solving the Implicit Euler ODE with Boundary Conditions

OK so far?Yes, that's correct. Sorry for the typo. And yes, everything is clear so far. Keep going!Great. So now, we can rewrite the two equations in matrix form as follows:$$\begin{pmatrix}1 & -\Delta x \\-\Delta x \frac{k}{\epsilon} & 1-\Delta x \frac{1}{\epsilon}\end{pmatrix}\begin{pmatrix}y_1^{n+1} \\y_2^{n+1}\end{pmatrix}=\begin{pmatrix}y_1^{n} \\y_2^{n}\end{pmatrix}$$This is a
  • #1
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##.
Code:
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
 
Physics news on Phys.org
  • #2
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.
 
  • #3
When you say you are using implicit Euler, does that mean you are solving this problem iteratively? Let's see your difference equations.
 
  • #4
DrClaude said:
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##?
Chestermiller said:
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
Code:
% 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)
 
  • #5
joshmccraney said:
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
Code:
% 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?"
 
  • #6
Chestermiller said:
So, you are solving this split boundary value problem using the "shooting method?"
What's a split bvp? And yep!
 
  • #7
You seem to be using a 2nd order central difference finite difference approximation when Euler is a first order difference.
 
  • #8
Chestermiller said:
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?
 
  • #9
joshmccraney said:
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?
 
  • #10
Chestermiller said:
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?
 
  • #11
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?
 
  • #12
Chestermiller said:
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 :)
 
  • #13
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?
 
Last edited:
  • #14
Chestermiller said:
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!
 
  • #15
joshmccraney said:
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.
 
  • #16
Chestermiller said:
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?
 
  • #17
joshmccraney said:
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?
 
  • #18
Chestermiller said:
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?
 
  • #19
joshmccraney said:
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.
 
  • #20
Chestermiller said:
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?
 
  • #21
joshmccraney said:
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.
 
  • Like
Likes member 428835
  • #22
Chestermiller said:
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?
 
  • #23
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.
 
  • #24
Thanks a lot!
 
  • #25
joshmccraney said:
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).
 

Related to Solving the Implicit Euler ODE with Boundary Conditions

1. What is the "Implicit Euler ODE BC" method?

The Implicit Euler ODE BC method is a numerical method used to solve ordinary differential equations (ODEs) with boundary conditions. It is an implicit method, meaning that it calculates the solution at the next time step using an equation that involves the solution at the current time step. This method is commonly used in scientific and engineering applications for its stability and accuracy.

2. How does the Implicit Euler ODE BC method differ from the Explicit Euler method?

The main difference between the Implicit Euler ODE BC method and the Explicit Euler method is the way they calculate the solution at the next time step. The Implicit Euler method uses an equation that involves the solution at the current time step, while the Explicit Euler method uses an equation that only involves the solution at the previous time step. This makes the Implicit Euler method more accurate and stable, but also more computationally expensive.

3. What are boundary conditions in the context of the Implicit Euler ODE BC method?

Boundary conditions are constraints that are imposed on the solution of an ODE at certain points or intervals. In the context of the Implicit Euler ODE BC method, boundary conditions are used to determine the values of the solution at the boundaries of the domain, which are typically specified as initial conditions or fixed values. These conditions are important for obtaining an accurate and meaningful solution to the ODE.

4. What are the advantages of using the Implicit Euler ODE BC method?

The Implicit Euler ODE BC method has several advantages over other numerical methods for solving ODEs. It is unconditionally stable, meaning that it can handle stiff equations and large time steps without the solution becoming unstable. It also has a high level of accuracy, making it suitable for solving complex and sensitive problems. Additionally, the Implicit Euler method can handle a variety of boundary conditions and is relatively easy to implement.

5. What are some common applications of the Implicit Euler ODE BC method?

The Implicit Euler ODE BC method has a wide range of applications in various fields of science and engineering. It is commonly used in physics, chemistry, biology, and engineering to model and simulate complex systems and phenomena. Some specific applications include population dynamics, chemical reactions, heat transfer, and fluid dynamics. The method is also used in computer graphics and animations to simulate motion and deformations of objects.

Similar threads

  • Calculus and Beyond Homework Help
Replies
6
Views
940
  • Calculus and Beyond Homework Help
Replies
10
Views
681
  • Calculus and Beyond Homework Help
Replies
12
Views
1K
  • Calculus and Beyond Homework Help
Replies
3
Views
649
  • Calculus and Beyond Homework Help
Replies
3
Views
481
  • Calculus and Beyond Homework Help
Replies
24
Views
1K
  • Calculus and Beyond Homework Help
Replies
2
Views
829
  • Calculus and Beyond Homework Help
Replies
18
Views
1K
  • Calculus and Beyond Homework Help
Replies
15
Views
1K
  • Calculus and Beyond Homework Help
Replies
2
Views
602
Back
Top