MATLAB 1D 2nd-Order nonlinear differential eqn using Matlab

[jam_27]

1. What is the physical problem?

2. Can you post your solution (a JPEG plot) to this particular problem?

3. Can you post your solution to the original problem as well as the precise form of the ODE you used and the boundary conditions?

4. Can I access via the net, the commercial FEM software you're using?


Hey Salty,

I hope this helps:
1) The physical problem is diffusion of electrons in a semiconductor laser.
2) Check attachment
3) Go through my first post , you will find a doc file for details & I have already given you the values of the constants.
4) The software’s are Ansys and FEMLAB. I have used both. Unfortunately they are not for free d/l. You can also use MATLAB's PDETOOL, if you have it.

Thanks for taking interest.

PS: Please do tell me the recipe that you are using. My recipe is given in the doc file. I think you are using a different one.

jam

 

[saltydog]

Thanks for that plot Jam. I'll try and work towards it.

Yea, your *.doc reference above. That's part of my problem I suppose. I don't have Microsoft Office so I couldn't ever view it. But that's Ok. I'll figure out how to do so one way or another.

Edit: Ok I downloaded a word viewer from the Microsoft site. I'll proceed . . .

The algorithm (recepie) that I use implements the method explained in the reference you gave above.

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

I believe I can obtain your result but need to get a better handle on high precision in Mathematica.

 

[saltydog]

Jam, I do believe I'm making progress with this. I noticed in your doc file that the first coefficient was negative. That changes things. This is the precise equation I solved:

y^{''}=-\frac{9.375}{3}\cdot 10^6+\frac{0.2}{3}y+\left(\frac{1.4}{3}\cdot 10^{-7}\right)y^2+\left(\frac{1}{3}\cdot10^{-9}\right)y-\frac{1}{3}\cdot 10^{-3};\quad y(-100)=y(100)=0

I've attached a plot of the results. As you can see it's not exactly like your solution. Is there a chance I'm not using the same BVP you used to generate that plot? Please let me know.

Thanks,
Salty

 

[jam_27]

Ok ya, I am sorry about the -ve sign. Its important and I missed it while explaining things to you. Sorry for that. And yes, I forgot to tell you one more thing. The first coefficient (with the -ve sign) on the RHS actually varies from x=-50 to x=50. That’s why I get the plot the way I get it. I am sure you will get the same plot now. By the way your solution looks perfect to me.

PS: Request: can you please explain in your own words the recipe you are using as per the following link. I didn't quite get the hang of it.
http://homepage.univie.ac.at/Franz.Vesely/cp0102/dx/node74.html

Thanks,
jam

 

[saltydog]

Ok ya, I am sorry about the -ve sign. Its important and I missed it while explaining things to you. Sorry for that. And yes, I forgot to tell you one more thing. The first coefficient (with the -ve sign) on the RHS actually varies from x=-50 to x=50. That’s why I get the plot the way I get it. I am sure you will get the same plot now. By the way your solution looks perfect to me.

PS: Request: can you please explain in your own words the recipe you are using as per the following link. I didn't quite get the hang of it.
http://homepage.univie.ac.at/Franz.Vesely/cp0102/dx/node74.html

Thanks,
jam

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&#039;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 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:<br /> <br /> \mathbf{A}\mathbf{Y}=-\mathbf{E}<br /> <br /> The matrix of the coefficients of the delta y terms is tridiagonal because each e&#039;i depends only on its nearest neighbors, the y+1, y, and y-1 variables. That&#039;s why the non-zero terms in the matrix go down diagonally. Right?<br /> <br /> 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&#039;s why iteration is necessary. Hope that helps. If not, which precise part of the algorithm is causing problems for you to understand?<br /> <br /> Also, in regards to:<br /> <br /> 1) The physical problem is diffusion of electrons in a semiconductor laser.<br /> <br /> 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.

 

[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}) &amp;= h(y_{i+1}+\Delta y_{i+1},y_i+\Delta y_i,y_{i-1}+\Delta y_{i-1}) \\ &amp;=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} \\ &amp;= e_i+\alpha_i \Delta y_{i-1}+ \beta_i \Delta y_{i}+ \gamma_i \Delta y_{i+1} \\ &amp;= 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&#039;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 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:<br /> <br /> \mathbf{A}\mathbf{Y}=-\mathbf{E}<br /> <br /> The matrix of the coefficients of the delta y terms is tridiagonal because each e&#039;i depends only on its nearest neighbors, the y+1, y, and y-1 variables. That&#039;s why the non-zero terms in the matrix go down diagonally. Right?<br /> <br /> 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&#039;s why iteration is necessary. Hope that helps. If not, which precise part of the algorithm is causing problems for you to understand?<br /> <br /> Also, in regards to:<br /> <br /> 1) The physical problem is diffusion of electrons in a semiconductor laser.<br /> <br /> 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.

<br /> <br /> <br /> Thanks for all the help.<br /> <br /> Experimentally we do it the following way.<br /> 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. <br /> Here we are talking of the simplest type of semiconductor laser, i.e. a broad area laser. <br /> <br /> We do the measurement below threshold (spontaneous emission) because above threshold (stimulated emission &amp; lasing) the situation is extremely complicated, which I am actually trying to model.<br /> <br /> This BVP is just a part of a bigger code which I am hope full will give me the results I am looking for.<br /> <br /> I hope I have been able to explain things properly.<br /> <br /> Now, the things you have explained about your recipe: I could already understand these stuff from the link:<br /> <a href="http://homepage.univie.ac.at/Franz.Vesely/cp0102/dx/node74.html" target="_blank" class="link link--external" rel="nofollow ugc noopener">http://homepage.univie.ac.at/Franz.Vesely/cp0102/dx/node74.html</a><br /> 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?<br /> g(M-1)=-a(M)/beta(M);<br /> h(M-1)=-e(M)/beta(M);<br /> g(i-1)=-alpha(i)/(beta(i)+gamma(i)*g(i));<br /> h(i-1)=(-e(i)-h(i))/(beta(i)+g(i));<br /> <br /> How does he get dely1=e1;<br /> <br /> May b the answers are trivial, but you see my math’s is very poor. I request you to kindly guide me.<br /> <br /> Thanks,<br /> <br /> 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:

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

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.