Why Does Increasing Final Run Time Cause Divergence in Finite Differencing?

[Chestermiller]

I do, since $h$ is even we evaluate from $[0,L]$ and multiply by $2$. Then we are left with the expression you wrote since $h_z(0,t)=0$. The last equality holds since volume is conserved. This is getting fun; what's the next plan? Obviously this qualifies your taylor expansion and post 24.

Yes, we are going to use Taylor series expansion. But, first I would like to go back to the transformed version of the problem:
$$\frac{\partial h^2}{\partial t}=Z\left(\frac{\partial h^2}{\partial Z}\right)\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{\partial ^2h^3}{\partial Z^2}\tag{1}$$
We note from the final equation of post #34 that volume will be conserved if, in the vicinity of Z=1, h is directly proportional to (Z-1). This also satisfies the boundary condition that h = 0 at Z = 1. So it makes sense to expand h in a Taylor series about Z = 1 as follows:
$$h=c_1(Z-1)+c_2(Z-1)^2+...\tag{2}$$
where
$$c_1(t)=\left(\frac{\partial h}{\partial Z}\right)_{Z=1}\tag{3}$$
$$c_2(t)=\frac{1}{2}\left(\frac{\partial^2 h}{\partial Z^2}\right)_{Z=1}\tag{4}$$
Retaining only the first two terms in Eqn. 2, what is h² and h³ (also retaining only the first two terms in these expressions)?

Chet

 

[member 428835]

So it makes sense to expand h in a Taylor series about Z = 1 as follows:
$$h=c_1(Z-1)+c_2(Z-1)^2+...\tag{2}$$
where
$$c_1(t)=\left(\frac{\partial h}{\partial Z}\right)_{Z=1}\tag{3}$$
$$c_2(t)=\frac{1}{2}\left(\frac{\partial^2 h}{\partial Z^2}\right)_{Z=1}\tag{4}$$
Retaining only the first two terms in Eqn. 2, what is h² and h³ (also retaining only the first two terms in these expressions)?
Chet

I assume you are looking for all terms before $\mathcal{O}(Z-1)^3$? If so I arrive at $$h^2 = c_1(t)(Z-1)^2 + \mathcal{O}(Z-1)^3 \\h^3 = 0 + \mathcal{O}(Z-1)^3$$
How does this look?

 

[Chestermiller]

I assume you are looking for all terms before $\mathcal{O}(Z-1)^3$? If so I arrive at $$h^2 = c_1(t)(Z-1)^2 + \mathcal{O}(Z-1)^3 \\h^3 = 0 + \mathcal{O}(Z-1)^3$$
How does this look?

I was thinking more like:

$$h^2=c_1^2(Z-1)^2+2c_1c_2(Z-1)^3$$
$$h^3=c_1^3(Z-1)^3+3c_1^2c_2(Z-1)^4$$

Next, please substitute these into the differential equation, and then write out the equations that are the coefficients of (Z-1) and (Z-1)².

Chet

 

[member 428835]

Ahh, now I se why you wanted the first two terms. We have $$2 c_1 c_1'(Z-1)^2 + 2(Z-1)^3(c_2c'_1+c_1 c_2') = Z\left[ 2 c_1^2(Z-1) + 6c_1 c_2 (Z-1)^2 \right]\frac{1}{L}\frac{dL}{dt} + \frac{2}{3L^2}\left[ 6c_1^3(Z-1)+36c_1^2c_2(Z-1)^2 \right] \implies\\
\mathcal{O}(Z-1): 0 = 3 L Z \frac{dL}{dt} + 2c_1\\
\mathcal{O}(Z-1)^2: c_1' L^2 = 3Zc_2 L\frac{dL}{dt} + 12c_1c_2 $$
Are you getting the same? I should say I have simplified and also canceled some $c_1$ terms since we are looking at a limit and can cancel.

 

[Chestermiller]

Yes. This expansion is accurate in the vicinity of Z = 1. If we now collect the coefficients of (Z-1), we get:
$$\frac{2Z c_1^2}{L}\frac{dL}{dt} + \frac{4c_1^3}{3L^2}=0$$
Rearranging this and evaluating it at Z = 1 (otherwise we could have written Z = 1+(Z-1)), we get
$$L\frac{dL}{dt}=-\frac{2}{3}c_1$$
So, if we know $c_1(t)$, we can integrate this to get L(t). Care to offer a finite difference approximation for $c_1$ in terms of h(1-ΔZ)?

Chet

 

[member 428835]

Care to offer a finite difference approximation for $c_1$ in terms of h(1-ΔZ)?
Chet

$$c_1(t) = \left( \frac{\partial h}{\partial Z} \right)_{Z=1} \approx \frac{h(t,1)-h(t,1-\Delta Z)}{\Delta Z}$$ yet $h(t,Z=1)=0$ thus we may write $$c_1(t) \approx \frac{h(t,1-\Delta Z)}{\Delta Z}$$
Is this what you were thinking?

 

[Chestermiller]

$$c_1(t) = \left( \frac{\partial h}{\partial Z} \right)_{Z=1} \approx \frac{h(t,1)-h(t,1-\Delta Z)}{\Delta Z}$$ yet $h(t,Z=1)=0$ thus we may write $$c_1(t) \approx \frac{h(t,1-\Delta Z)}{\Delta Z}$$
Is this what you were thinking?

Yes, except for a missing minus sign. Now you can substitute this into the equation for dL/dt to obtain ?

Chet

 

[member 428835]

Yes, except for a missing minus sign. Now you can substitute this into the equation for dL/dt to obtain ?
Chet

We have $$ \int_{L_0}^L L' \, dL' \approx \frac{2}{3} \int_0^t \frac{h(t',1-\Delta Z)}{\Delta Z} \, dt' \implies \\
L = \sqrt{L_0^2 + \frac{4}{3} \int_0^t \frac{h(t',1-\Delta Z)}{\Delta Z} \, dt'}$$ where prime notation denotes dummy variables and $L_0$ is the initial length. I dropped the minus sign after the root since we can work with positive $Z$ domain. How to solve the integral though?

 

[Chestermiller]

$$L^2(t+Δt)=L^2(t)+\frac{4}{3}\frac{h(t,1-ΔZ)}{\Delta Z}Δt$$
In the main differential equation, you need L², and you need $\frac{1}{L}\frac{dL}{dt}$. But,
$$\frac{1}{L}\frac{dL}{dt}=\frac{L}{L^2}\frac{dL}{dt}=\frac{1}{2L^2}\frac{dL^2}{dt}$$
This is approximated numerically by:
$$\frac{1}{2L^2}\frac{dL^2}{dt}≈\frac{\frac{2}{3}\frac{h(t,1-ΔZ)}{\Delta Z}}{L^2(t)+\frac{2}{3}\frac{h(t,1-ΔZ)}{\Delta Z}Δt}$$
Note that this is evaluated half way between t and t + Δt (for better numerical accuracy). In the main differential equation, you also evaluate the L² at t+Δt/2:
$$L^2(t+Δt/2)=L^2(t)+\frac{2}{3}\frac{h(t,1-ΔZ)}{\Delta Z}Δt$$
So, at each time step, you first evaluate the new values of the two expressions involving L in the main differential equation and then use these same values at each grid point to evaluate the new h^2's.

Chet

 

[member 428835]

$$L^2(t+Δt)=L^2(t)+\frac{4}{3}h(t,1-ΔZ)Δt$$

I'm trying to follow your work, and it looks like you're doing this: $$ \int_{L(t)}^{L(t+\Delta t)} L dL' = \frac{2}{3}\int_{t}^{t+\Delta t}\frac{h(t',1-\Delta Z)}{\Delta Z} dt' $$ Evidently $$\int_{t}^{t+\Delta t}\frac{h(t',1-\Delta Z)}{\Delta Z} dt' = h(t',1-\Delta Z) \Delta t$$ Could you explain this? I don't see how the $\Delta Z$ is gone.

This is approximated numerically by:
$$\frac{1}{2L^2}\frac{dL^2}{dt}≈\frac{\frac{2}{3}h(t,1-ΔZ)Δt}{L^2(t)+\frac{2}{3}h(t,1-ΔZ)Δt}$$

I'm also having trouble here, if you can help? I understand that $L^2$ in the denominator is $L^2(t)+\frac{2}{3}h(t,1-ΔZ)Δt$ from above (although shouldn't it be $L^2(t)+\frac{4}{3}h(t,1-ΔZ)Δt$?) If so, this means $$\frac{1}{2}\frac{d L^2}{dt} = \frac{2}{3} h(t,1-\Delta Z)\Delta t$$. Could you also explain this please?

Thanks so much! You're great!

 

[Chestermiller]

I'm trying to follow your work, and it looks like you're doing this: $$ \int_{L(t)}^{L(t+\Delta t)} L dL' = \frac{2}{3}\int_{t}^{t+\Delta t}\frac{h(t',1-\Delta Z)}{\Delta Z} dt' $$ Evidently $$\int_{t}^{t+\Delta t}\frac{h(t',1-\Delta Z)}{\Delta Z} dt' = h(t',1-\Delta Z) \Delta t$$ Could you explain this? I don't see how the $\Delta Z$ is gone.

Oops. I left out the delta z. I'm going back and correcting all the equations. Thanks for spotting this.

I'm also having trouble here, if you can help? I understand that $L^2$ in the denominator is $L^2(t)+\frac{2}{3}h(t,1-ΔZ)Δt$ from above (although shouldn't it be $L^2(t)+\frac{4}{3}h(t,1-ΔZ)Δt$?)

No. It's being evaluated at the time half-way across the interval to provide better accuracy in the time integration of the main differential equation.

If so, this means $$\frac{1}{2}\frac{d L^2}{dt} = \frac{2}{3} h(t,1-\Delta Z)\Delta t$$. Could you also explain this please?

Ooooops. Another mistake. Going back to correct that too. Please give me about 10 minutes to correct it, and then let me know if it's OK.

Chet

 

[member 428835]

I think I follow you. Let me reiterate it and you let me know if I've made a mistake. Begin with

$$\int_{L(t)}^{L(t+\Delta t)}L' dL' =\frac{2}{3} \int_t^{t+\Delta t}\frac{h(t',1-\Delta Z)}{\Delta Z} dt' \implies\\
L^2(t+\Delta t) = L^2(t) + \frac{4}{3} \frac{h(t',1-\Delta Z)}{\Delta Z} \Delta t $$ Thus we have an expression for $L^2(t+\Delta t)$. The integral was approximated using a left endpoint rule. Next we must find $dL/dt$. Begin where we left off

$$ L^2(t+\Delta t) = L^2(t) + \frac{4}{3} \frac{h(t',1-\Delta Z)}{\Delta Z} \Delta t \implies\\
\frac{L^2(t+\Delta t) - L^2(t)}{\Delta t} = \frac{4}{3} \frac{h(t',1-\Delta Z)}{\Delta Z}$$ Notice the left side term is an approximation of $dL/dt$. Plugging in the above into
$$\frac{1}{L}\frac{dL}{dt}=\frac{L}{L^2}\frac{dL}{dt}=\frac{1}{2L^2}\frac{dL^2}{dt}$$
Is approximated numerically by:
$$\left.\frac{1}{2L^2}\frac{dL^2}{dt}\right|_{t=t+\Delta t} = \frac{\frac{2}{3} \frac{h(t',1-\Delta Z)}{\Delta Z}}{L^2(t) + \frac{4}{3} \frac{h(t',1-\Delta Z)}{\Delta Z} \Delta t }$$

I'm missing the $t+\Delta t/2$ concept, though.

 

[Chestermiller]

I think I follow you. Let me reiterate it and you let me know if I've made a mistake. Begin with

$$\int_{L(t)}^{L(t+\Delta t)}L' dL' =\frac{2}{3} \int_t^{t+\Delta t}\frac{h(t',1-\Delta Z)}{\Delta Z} dt' \implies\\
L^2(t+\Delta t) = L^2(t) + \frac{4}{3} \frac{h(t',1-\Delta Z)}{\Delta Z} \Delta t $$

I would replace the t primes in the finite difference equation with t. That goes for all the other finite difference equations below as well.

Thus we have an expression for $L^2(t+\Delta t)$. The integral was approximated using a left endpoint rule. Next we must find $dL/dt$. Begin where we left off

$$ L^2(t+\Delta t) = L^2(t) + \frac{4}{3} \frac{h(t',1-\Delta Z)}{\Delta Z} \Delta t \implies\\
\frac{L^2(t+\Delta t) - L^2(t)}{\Delta t} = \frac{4}{3} \frac{h(t',1-\Delta Z)}{\Delta Z}$$ Notice the left side term is an approximation of $dL/dt$. Plugging in the above into
$$\frac{1}{L}\frac{dL}{dt}=\frac{L}{L^2}\frac{dL}{dt}=\frac{1}{2L^2}\frac{dL^2}{dt}$$
Is approximated numerically by:
$$\left.\frac{1}{2L^2}\frac{dL^2}{dt}\right|_{t=t+\Delta t} = \frac{\frac{2}{3} \frac{h(t',1-\Delta Z)}{\Delta Z}}{L^2(t) + \frac{4}{3} \frac{h(t',1-\Delta Z)}{\Delta Z} \Delta t }$$

I'm missing the $t+\Delta t/2$ concept, though.

For extra accuracy, I'm trying to evaluate the entire term at the average over the time interval between t and t + Δt, rather than at either end of the interval. It really doesn't matter much whether, in the denominator, you use L²(t), L²(t+Δt/2), or L²(t+Δt). L²(t+Δt/2) just gives a little more accuracy at no additional computational cost. If you want to use the value at t + Δt, I have no major problem with that. I just feel that, if the added accuracy is there for the taking, why not take it?

Chet

 

[member 428835]

I would replace the t primes in the finite difference equation with t. That goes for all the other finite difference equations below as well.

Ooops, yea, I didn't see this. My fault.

For extra accuracy, I'm trying to evaluate the entire term at the average over the time interval between t and t + Δt, rather than at either end of the interval.
Chet

So it seems what changes is the 4 in the denominator turns into a 2, right? Does this stem from $$\frac{\partial h}{\partial t} \approx \frac{h(t+\Delta t/2)-h(t)}{\Delta t/2}$$

So then we are solving $$\left(\frac{\partial y}{\partial t}\right)_Z=Z\left(\frac{\partial y}{\partial Z}\right)_t\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{\partial ^2(y^{3/2})}{\partial Z^2}$$

where we just found an expression for the $L$ terms. Boundary conditions are $y(t,Z=\pm 1)=0$, and $y'(t,Z=0)=0$. Am I missing any BC's and how would we take a finite difference of the non-linear term? We can't just write it as $(y_i^j)^{3/2}$ can we? The chain rule should be considered, right? Or should we not make the $y$ transform and stay with $h$ using instead the $Z=z/L(t)$?

 

[Chestermiller]

Ooops, yea, I didn't see this. My fault.
So it seems what changes is the 4 in the denominator turns into a 2, right? Does this stem from $$\frac{\partial h}{\partial t} \approx \frac{h(t+\Delta t/2)-h(t)}{\Delta t/2}$$

No. We are still doing
$$\frac{\partial y}{\partial t} \approx \frac{y(t+\Delta t)-y(t)}{\Delta t}$$
But you get extra accuracy if you evaluate the right hand side of the equation at t + Δt/2.

So then we are solving $$\left(\frac{\partial y}{\partial t}\right)_Z=Z\left(\frac{\partial y}{\partial Z}\right)_t\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{\partial ^2(y^{3/2})}{\partial Z^2}$$

where we just found an expression for the $L$ terms. Boundary conditions are $y(t,Z=\pm 1)=0$, and $y'(t,Z=0)=0$. Am I missing any BC's

Yes. I'm hoping that you are only solving the equation over the domain from 0 to 1, and not from -1 to + 1. And I'm hoping that, at y = 0, you are solving the differential equation for y(0) using a reflective boundary condition, rather than just setting y(0) equal to y(Δz)

and how would we take a finite difference of the non-linear term? We can't just write it as $(y_i^j)^{3/2}$ can we?

Yes. This is what you should do.

The chain rule should be considered, right? Or should we not make the $y$ transform and stay with $h$ using instead the $Z=z/L(t)$?

No and no.

Chet

 

[member 428835]

I'm hoping that you are only solving the equation over the domain from 0 to 1, and not from -1 to + 1.

Absolutely!

And I'm hoping that, at y = 0, you are solving the differential equation for y(0) using a reflective boundary condition, rather than just setting y(0) equal to y(Δz)

I'm a little unsure of what you are saying here. Can you specify please?

 

[Chestermiller]

Absolutely!
I'm a little unsure of what you are saying here. Can you specify please?

$$\left(\frac{\partial y}{\partial t}\right)_Z=Z\left(\frac{\partial y}{\partial Z}\right)_t\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{\partial ^2(y^{3/2})}{\partial Z^2}$$
At z = 0, $\frac{\partial y}{\partial Z}=0$, then $y^{3/2}(-ΔZ)=y^{3/2}(+ΔZ)$. This means that:
$$y^{3/2}(-ΔZ)-2y^{3/2}(0)+y^{3/2}(+Δz)=2[y^{3/2}(+Δz)-y^{3/2}(0)]$$
So, at Z = 0,
$$\frac{y(0,t+Δt)-y(0,t)}{Δt}=\frac{4}{3L^2}\frac{[y^{3/2}(t,+Δz)-y^{3/2}(t,0)]}{(ΔZ)^2}$$

This is the finite difference equation that should be used at Z = 0.

Chet

 

[member 428835]

I'm with you! Thanks for clearing that up! Ok, so since our thoughts have been all over the place (mainly due to my questions) I've decided to type up a latex document. I'll post the pdf when I finish it for both your review and others who may want to see the answer but don't want to look through all the conversation. Thanks a ton for your help. I'll probably be writing on this thread before the pdf is finished, as I want to ensure the FD is correct.

 

[Chestermiller]

Please note also, that, as things stand now, you are still using Forward Euler to integrate with respect to time. If stability is an issue, you can control it by taking the time step small enough. But please note that, in terms of stability, Forward Euler is the worst scheme that one can use. So, if the time step is too large, the scheme will go unstable. Stability problems can typically be removed by going to so-called implicit schemes, such as Backward Euler, in which the right hand side of the equation is evaluated at t + Δt, rather than t. But this requires the solution of a set of coupled algebraic equations at each time step. Hopefully, for your problem, Forward Euler will be adequate. But, if not, I can help.

Chet

 

[member 428835]

Chet, I have a quick question about taking a finite difference over non-linear terms. Notice the following $$ \frac{\partial (h^2)}{\partial t} = 2h\frac{\partial h}{\partial t}$$ Applying a finite difference to each yields $$\frac{\partial (h^2)}{\partial t} = \frac{h^2(t+\Delta t) - h^2(t)}{\Delta t}$$ but also $$ 2h\frac{\partial h}{\partial t} = 2 h(t) \frac{h(t+\Delta t)-h(t)}{\Delta t} $$ Clearly the finite differences are not equal. Can you explain what is happening?

 

[member 428835]

Attached is the finite difference technique along with how we arrived. Please review it, especially the end, and let me know what you think. The conclusion doesn't look right. EDIT: the limits should be from the left, but you get the idea.

 

[Chestermiller]

Chet, I have a quick question about taking a finite difference over non-linear terms. Notice the following $$ \frac{\partial (h^2)}{\partial t} = 2h\frac{\partial h}{\partial t}$$ Applying a finite difference to each yields $$\frac{\partial (h^2)}{\partial t} = \frac{h^2(t+\Delta t) - h^2(t)}{\Delta t}$$ but also $$ 2h\frac{\partial h}{\partial t} = 2 h(t) \frac{h(t+\Delta t)-h(t)}{\Delta t} $$ Clearly the finite differences are not equal. Can you explain what is happening?

These are two slightly different first order finite difference approximations to the exact same quantity, and, as such, neither one of them is favored over the other. In each case, the discretization error decreases linearly with Δt when Δt becomes small. To see that, expand h(t+Δt) in a Taylor series about t, and substitute the resulting series into each of the two finite difference approximations.

Chet

 

[Chestermiller]

Attached is the finite difference technique along with how we arrived. Please review it, especially the end, and let me know what you think. The conclusion doesn't look right. EDIT: the limits should be from the left, but you get the idea.

Corrections:

Eqn. 16 should be the derivative of h³, not h.

In Eqn. 25, there should be a 2 in the numerator, not a 4. I would also prefer to see a 4 in the denominator, but I can't strongly argue against the 4.

Eqn. 26 is incorrect. The second derivative of h³ with respect to Z has been evaluated incorrectly.

I didn't check Eqns. 28 and 29.

Josh,

I've had a huge amount of practical experience with numerical analysis, particularly on problems like this one. And that practical experience involved many situations in which a lot of money was riding on the outcome of the correct numerical solution of the equations. Based on this experience, I am very strongly urging you to solve this problem numerically in terms of y rather than in terms of h, for the following reasons:

1. There are many fewer mathematical operations involved in solving it this way
2. There is much less of a chance of making a coding error
3. The finite difference equations that result in this approach automatically conserve volume, unlike the scheme involving h, in which volume is conserved only in the limit of vanishing $\Delta Z$

If I were coding this problem, starting from the beginning of each time step (so that y is known at all grid points at time t), I would carry out the following sequence of steps:

1. Calculate y^3/2 at all the grid points
2. Calculate √y at the grid point adjacent to the right boundary. This is h at that location.
3. Calculate L²(t+Δt) and (1/L)dL/dt
4. Evaluate the right hand side of the main differential equation at each grid point
5. Execute a time step at each grid point

Chet

 

[member 428835]

These are two slightly different first order finite difference approximations to the exact same quantity, and, as such, neither one of them is favored over the other. In each case, the discretization error decreases linearly with Δt when Δt becomes small. To see that, expand h(t+Δt) in a Taylor series about t, and substitute the resulting series into each of the two finite difference approximations.
Chet

Ok, this is cool! I just did the expansion and it works. Awesome!

 

[member 428835]

Josh,

I've had a huge amount of practical experience with numerical analysis, particularly on problems like this one. And that practical experience involved many situations in which a lot of money was riding on the outcome of the correct numerical solution of the equations. Based on this experience, I am very strongly urging you to solve this problem numerically in terms of y rather than in terms of h
Chet

I'll definitely do it this way now that I see (from your previous post) that the expansion is the same! This was not immediately obvious to me, hence why I never introduced the $y$ transform. Thanks so much fro your help, Chet! Stay tuned, as I'll fix the corrections, introduce the transform, and also post a flow chart for the code (although you've seem to have already done this).

You're awesome!

 

[member 428835]

Ok, I think I corrected all the errors you pointed out. How does this attached version look?

As an aside, do we really need (12) through (16)? They seem to not really be used.

 

[Chestermiller]

You're awesome!

Not really. Just lots of experience.

 

[Chestermiller]

Ok, I think I corrected all the errors you pointed out. How does this attached version look?

As an aside, do we really need (12) through (16)? They seem to not really be used.

Yes. I would lose 12 - 16.

There's a typo in Eqn. 29: the y should be under the radical in the numerator.

Please show the derivation of Eqns. 30 and 34-36. I'm not able to follow in detail. I just want to make sure that it was done correctly. For example, I think there should be a 1/2 in Eqn. 30 instead of the factor of 2; and, in Eqn. 34, I feel like there should be a 4 in the denominator. Also, in Eqn. 36, there should be a ζ in the denominator rather than a ξ.

Chet

 

[member 428835]

How does this look?

(23) is evaluated at $t$ and yet requires $y$ to be evaluated at $t$ as well. How would you deal with this?

I also got rid of zeta and xi.

 

[Chestermiller]

How does this look?

I also got rid of zeta and xi.

I found a couple of errors in the new version:
Eqns. 17,19, 22,23,33,34,and 35 should have a 4 instead of a 2.

In Eqn. 32, you are missing an additional factor of 2 in the denominator of the first term in brackets. This is because ∂y/∂Z is evaluated numerically using a central difference over 2 grid intervals.

(23) is evaluated at $t$ and yet requires $y$ to be evaluated at $t$ as well. How would you deal with this?

In Eqns. 33-34, you need to replace the j's with (j+1) and you need to replace the (j-1)'s with j. And in Eqn. 35, you need to replace the (j-1) with j (so you have j on both sides). As far as the L_j² in the denominator of Eqn. 32 is concerned, you have already calculated that from the previous time step (Or, to the same degree of accuracy, you could also use L_j+1² or the average of the two).

 

[member 428835]

Thanks for pointing these errors out. I agree with all of your corrections, and believe this is now error free. Thanks for all the help!

How would I go about an error analysis though? It's non-linear so I don't think a von Neumann approach is valid since separation of variables doesn't work.

I'll start coding this too (in MatLab) and post a few pictures of what I'm getting in the .pdf for you're consideration.

 

[Chestermiller]

Thanks for pointing these errors out. I agree with all of your corrections, and believe this is now error free. Thanks for all the help!

How would I go about an error analysis though? It's non-linear so I don't think a von Neumann approach is valid since separation of variables doesn't work.

I'll start coding this too (in MatLab) and post a few pictures of what I'm getting in the .pdf for you're consideration.

Are you talking about error analysis, or are you talking about stability?

Chet

 

[member 428835]

Are you talking about error analysis, or are you talking about stability?

Chet

Sorry, I'm talking about stability.

 

[Chestermiller]

Sorry, I'm talking about stability.

I don't know of any methods of analyzing the stability of time-dependent schemes in non-linear problems, but, experience has shown that stability is guaranteed if you go to fully implicit (e.g., Backward Euler). It is also possible to increase stability by going partially implicit (i.e., more implicit) by handling certain quantities implicitly, while treating others explicitly.

Solving the problem implicitly means that you have to solve a set of non-linear equations at each time step. The easiest method of solving a problem like this implicitly is to use the so-called "method of lines" in conjunction with a stiff ODE software package based on the C. W. Gear method. Such a package transparently solves the set on non-linear equations at each time step. The method of lines changes the non-linear PDE into a non-linear set of coupled ODEs, by finite differencing the spatial derivatives, but not the time derivative. So there is an ODE for each grid point. The stiff package automatically integrates these stiff ODEs.

A partially implicit method can be developed for your problem that requires a tri-diagonal set of linear equations to be solved at each time step. This can be done using a tri-diagonal matrix inversion routine. However, there will be some quantities in the difference equations that will have to be handled explicitly in this approach. Still, in my judgement, it would improve stability very substantially (and allow much larger time steps).

The first thing to do is to see whether the existing scheme gives you stability when small (but acceptable) time steps are used.

Chet

 

[member 428835]

Thanks a ton! So I just finished writing the code. It's at the end of the attached pdf. I put plots in the pdf for a time value of 14 (line 8 in the code) and it ran well. However, if I set $tf = 15$ I get some errors, specifically imaginary parts. My guess is from line 37/38, when the square root is introduced. I should say that if I change line 9 to $N = (tf/10)*M^2;$ I can run for $tf$ values of 100 and more.

What is your take on this technique? Is this acceptable or should I go for an implicit scheme?

Also, for your curiosity, the percent change in volume was less than 1.5% for the code and pictures I posted. Changing line 9 as suggested and running for $tf=100$ produces a volume change of about 0.003%.