1D 2nd-Order nonlinear differential eqn using Matlab

[jam_27]

Good Jam. I do believe you've taught me something.

Jam, the method relies on the Taylor expansion of the error terms. Above in post #25, I showed how each error term for each value of y calculated depends on three variables:

$$e_i=F(y_{i+1},y_i,y_{i-1})$$

We wish to calculate the value of this function (for each $e_i$)
for small changes in the variables as per posting #25 above:

$$\begin{align*}e_i(\mathbf{Y}+\Delta\mathbf{Y}) &= h(y_{i+1}+\Delta y_{i+1},y_i+\Delta y_i,y_{i-1}+\Delta y_{i-1}) \\ &=e_i+\frac{\partial e_i}{\partial y_{i+1}} \Delta y_{i+1}+\frac{\partial e_i}{\partial y_{i}} \Delta y_{i}+\frac{\partial e_i}{\partial y_{i-1}} \Delta y_{i-1} \\ &= e_i+\alpha_i \Delta y_{i-1}+ \beta_i \Delta y_{i}+ \gamma_i \Delta y_{i+1} \\ &= e_i+\Delta y_{i-1}+\left(-2-(\Delta x)^2[-1+3(y_i)^2]\right)\Delta y_i+\Delta y_{i+1} \\\end{align}$$

We treat the $\Delta y_i's$ as unknowns and set up the n equations so that each is set to zero (that will make each error term zero). Since each term includes an [itex]e_i[/tex] term as per the taylor series of any function, we put those e terms on the right hand side and end up with the matrix equation:

$$\mathbf{A}\mathbf{Y}=-\mathbf{E}$$

The matrix of the coefficients of the delta y terms is tridiagonal because each e'i depends only on its nearest neighbors, the y+1, y, and y-1 variables. That's why the non-zero terms in the matrix go down diagonally. Right?

Solving this tridiagonal matrix equation gives us these values which we then add to the current values of y. Since there is a slight error associated with any taylor expansion, the results are not perfect. That's why iteration is necessary. Hope that helps. If not, which precise part of the algorithm is causing problems for you to understand?

Also, in regards to:

1) The physical problem is diffusion of electrons in a semiconductor laser.

Can you explain how you experimentally verified that the diffusion follows the solution to the BVP? For example, could the output be monitored by an O-scope or some other device which reflects the square-wave phenomenon? Just curious.

Thanks for all the help.

Experimentally we do it the following way.
We open a window on the n-side (substrate) of the semiconductor laser photolithographic ally and measure the spontaneous emission profile of the light, which theoretically should follow the carrier diffusion in the device. We expect a square current density profile (first term on the RHS (-ve)) and in accordance to that we expect almost square carrier diffusion in the device.
Here we are talking of the simplest type of semiconductor laser, i.e. a broad area laser.

We do the measurement below threshold (spontaneous emission) because above threshold (stimulated emission & lasing) the situation is extremely complicated, which I am actually trying to model.

This BVP is just a part of a bigger code which I am hope full will give me the results I am looking for.

I hope I have been able to explain things properly.

Now, the things you have explained about your recipe: I could already understand these stuff from the link:
http://homepage.univie.ac.at/Franz.Vesely/cp0102/dx/node74.html
What however what were unclear to me are the following things: why are these definitions needed. What purpose does these serve. Why is the author defining this stuff?
g(M-1)=-a(M)/beta(M);
h(M-1)=-e(M)/beta(M);
g(i-1)=-alpha(i)/(beta(i)+gamma(i)*g(i));
h(i-1)=(-e(i)-h(i))/(beta(i)+g(i));

How does he get dely1=e1;

May b the answers are trivial, but you see my math’s is very poor. I request you to kindly guide me.

Thanks,

Jam

 

[saltydog]

Thanks for all the help.

Experimentally we do it the following way.
We open a window on the n-side (substrate) of the semiconductor laser photolithographic ally and measure the spontaneous emission profile of the light, which theoretically should follow the carrier diffusion in the device. We expect a square current density profile (first term on the RHS (-ve)) and in accordance to that we expect almost square carrier diffusion in the device.
Here we are talking of the simplest type of semiconductor laser, i.e. a broad area laser.

We do the measurement below threshold (spontaneous emission) because above threshold (stimulated emission & lasing) the situation is extremely complicated, which I am actually trying to model.

Alright, it's over my head. Still interesting though.

What however what were unclear to me are the following things: why are these definitions needed. What purpose does these serve. Why is the author defining this stuff?
g(M-1)=-a(M)/beta(M);
h(M-1)=-e(M)/beta(M);
g(i-1)=-alpha(i)/(beta(i)+gamma(i)*g(i));
h(i-1)=(-e(i)-h(i))/(beta(i)+g(i));

How does he get dely1=e1;
Jam

Oh, now I see what you mean: That's the part where he actually solves the matrix equation. Well, as you know, I just kind of skidded around that part in my code by using Mathematica's built-in "TridiagonalSolve" function. I suspect the code you're having problems with above is the implementation of solving a "diagonally dominant" matrix equation via Gauss-Seidel iteration or Successive Overrelaxation Iteration. Ohhhh. Now I see why it's called "Relaxation Method" Suppose he's using SOR then. Might take me a little while of reviewing to figure out how that's being implemented above. I tell you what, the best way to understand this is to for the moment put away this problem and just work through a simple say 5x5 tridiagonal matrix equation using SOR. Then I bet we could understand what he's doing above.

 

[saltydog]

Tell you what Jam, I ain't proud. Let's solve this one using SOR. I bet the process will reflect what the web reference is doing:

$$\left( \begin{array}{cccc}4 & 5 & 0 & 0 \\ 7 & 8 & 9 & 0 \\ 0 & 1 & 2 & 3 \\ 0 & 0 & 11 & 12 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \\ x_4 \end{array} \right)= \left( \begin{array}{c} 2 \\ 3 \\ 4 \\ 5 \end{array} \right)$$

You're not offended I'm suggesting we look at such a baby equation are you? It's like they say: "sometimes you have to take two steps backwards to move forwards".

 

[jam_27]

I am glad that you are there to help me out. The only problem is that I would be busy for the next 2 days, i will get qualtiy time to look at only on Thursday. I will try to look at it today however.Lets see...

Thanks,

Jam

 

[jam_27]

And yes, I will have to look at SOR too. By the way this particular Matrix equation can be solved very easily (a 2 line code) in MATLAB. However let me look at SOR.I will come back to you.
Jam

 

[saltydog]

Jam, for the record, that reference you gave regarding the solution of a tridiagonal matrix equation is NOT using any of the three iterative methods (Jacobi, Gauss-Seidel, and SOR) but rather a direct substitution method. This is Ok for N not too large as it requires a lot more memory than the iterative techniques. I can understand why the tridiagonal matrix associated with numerical solutions to PDEs are solved iteratively since we're working on a 2-D (or higher dim) grid to solve them and thus have much more points than with an ODE.

Further up the pages on that website describes the recursive method in a little more detail:

http://homepage.univie.ac.at/Franz.Vesely/cp0102/dx/node23.html

Although I understand what's going on, it's "encapsulating" the details nicely in terms of those g and h variables which make it (for me anyway) tough to follow the notation.

 

[jam_27]

Thanks for the link , Salty...I will look into it.
Jam

 

[oasiszhu]

hi,saltydog. can you introduce some numerical method for dealing with the nonlinear partial differential equation with an initial value.