How can I run this iteration?

I want to calculate the pressure in a soil column for different depths i and different times t. My solution calculates the pressure u at time t + delta t for each node i sequentially, which means, I am looking for u(i,t+delta t). Initial pressures u(i,t) are given.

However, the solution of u(i,t+delta t) depends on u(i+1,t+delta t), which has not yet been calculated, so that the equation needs iterations.

All other variables (Cv, h) are known as well and do not change.

I don't know how to implement this in Matlab. Can anyone help me? It would help me so much!

Thanks a lot,

Juditha

Answer 1

[Bjorn Gustavsson]

This look very similar to the Crank-Nicolson method of solving a partial differential equation. In the way you've phrased the equation you have your unknowns (u at the next time-step) on both sides. What you need to do is to move all unknown to the left-hand side. This will give you a left-hand-side matrix M_l that is multiplied with the unknowns, and a right-hand-side with a couple of terms:

Maybe you will have no , that's not critical. This system of equations you can solve for the next time-step:

or in pseudo-matlab code:

For the matrices you have to make space for your boundary conditions. This will take you "a day" to implement, start with pen and paper and write the equations for one component of u then the two neighbouring, implement it in matlab, run the problem for one time-step for a very small problem (10 altitudes), in debug-mode so that you can look at the matrices and the solution. Also have a look at the Crank–Nicolson method. If you so this make heavy use of built-in functions like diag to create the differential-operation matrices. The C-N method is stable for you type of equation, however there are numerical conditions on the time-step. So when integrating the pde this way you also have to read-up on the Courant-Friedrich-Levy condition that tells you how long time-steps to take for a given spatial resolutino.

You can also have a look at matlab's PDE-solver pdepe and try to rephrase the problem to use that function.

HTH

[Juditha Schmidt]

Thank you for this quick and detailed advice! I will have a look at it and tell you if it works.

[Juditha Schmidt]

I found a solution for my problem! Thanks a lot! That was great!

[Bjorn Gustavsson]

Good! Happy to have helped. Did you go for a C-N type solution?

[Juditha Schmidt]

Yes! The C-N method works perfectly in my case.