Coupled PDE System - Numerical Solution

[FrankST]

All,

As part of my research I came up with a boundary value problem where I need to solve the following system of coupled PDE:

1- a1 * f,xx + a2 * f,yy + a3 * g,xx + a4 * g,yy - a5 * f - a6 * g = 0

2- b1 * f,xx + b2 * f,yy + b3 * g,xx + b4 * g,yy - b5 * f - b6 * g = 0


Where, ai's and bi's are all positive constants.

To solve this system numerically, I discretize it using finite difference scheme. In Eq. 1, I move all the terms associated with g(x,y) to the RHS and I solve Eq.1 for f(x,y).

Then, in Eq.2 I move all terms associated with f(x,y) (that is updated from the previous step) to RHS and I solve Eq.2 for g(x,y) and I continue this process until f and g converge to a unique solution.

I am wondering if there exists a more efficient way than my method to solve this system numerically and using finite difference scheme.

Your comments are highly appreciated.


Frank

 

[meldraft]

Express it as an optimization problem and minimize it to zero.

If your two expressions are ex1 and ex2, then:

$$N(x,y)=ex1^2+ex2^2$$

You should be able to minimize N(x,y) to zero quite easily by using the routines built into Matlab or similar software. In fact, even steepest descent with a randomized starting point could do the trick. Your guide into solving the problem is the value of f that you find. Since the system is supposed to have a solution it should be something like N(x,y)=1.e-8. If you get e.g. N(x,y)=2, then you have hit a local minimum and you need to start from a different starting point.

Hope it helps ;)Edit: replaced f with N

 

[FrankST]

Thanks for your answer.

But I am not sure if I understood your method. I can't see any g(x,y) in your solution.

I have two unknown functions called f(x,y) and g(x,y). I need to find these two functions by solving the system of PDE I mentioned earlier.

My question was that if there is a more efficient way than the iterative method that I used to solve this system.

 

[meldraft]

Oh sorry I accidentally just used the same letter as in your problem. You need to define a new function, let's call it N(x,y) (the one I originally called f(x,y)). I will edit it in my previous post.

Thanks for your answer.
I have two unknown functions called f(x,y) and g(x,y). I need to find these two functions by solving the system of PDE I mentioned earlier.

This is not solving the system though. Solving the system would mean finding the x,y that satisfy both equations. Moreover, solving the system numerically would give you numbers, not expressions, so I'm now a little confused :uhh:

The solving speed will depend on the form of the f,g functions as well as how you write the code, and of course which optimization algorithm you use and the sensitivity of the problem to the starting point. Solving a somewhat more complex 5x5 system than this in MATLAB using the optimization algorithms takes considerably less than a second on a 2.8 i5 CPU.

 

[FrankST]

By solving a system I do mean finding f(x,y) and g(x,y) that satisfy the PDE system I showed.

I have a rectangular domains with thousands of nodes (discretized domain). I need to find f and g at each node. Therefore, I am dealing with a system of equations (after discretization) consisting of thousands of equations.

 

[meldraft]

Thank you, this makes it much clearer. If you are solving it over a domain then it's a different story, I thought you had a symbolic expression and you wanted the solution at one point.

You could still replace your two equations with one ( N(x,y) ), and optimize over the domain until it converges. Depending on the problem, it might converge pretty fast (a Newton method may be the way to go in this case).

You'll probably avoid a lot of pain if you use a PDE solver though.

If you already know the stuff below, just ignore it ;)

If your geometry is not too complex, you can try using finite elements. It will probably require a much coarser mesh for the same level of accuracy. There are solvers freely available that don't require too much work to set up. One set of routines I regularly use is the routines in MATLAB's PDE toolbox, or, for more custom equations like yours, OpenFoam.

If you absolutely have to do it with finite differences, I'm afraid I can't help you too much as I only used these schemes while studying. What I do know is that there are 3 things that are of importance if it's taking too long:

1. The finite difference formula you are using for the level of accuracy you need
2. If your system has a RAM bottleneck or
3. a CPU bottleneck

If RAM usage peaks, check matrix memory allocation/deallocation. If you don't use all your CPU, try sending part of the calculations to a different thread/CPU or GPU.

 

[FrankST]

Thanks a lot for your valuable information.

I tried an analytical solution first but due to geometry of the domain (there are some holes inside the rectangular domain where I have Dirichlet BCs) I decided to use a numerical method and FDM is my current option.

In the iterative method that I am using it takes sometimes around 1000 iterations to converge to a solution (I set the tolerance as 1e-8). That's why I asked if there is a more efficient way to solve this system.

However, your hint on the FD formula and level of accuracy can be a good start for me.

Thanks for your time.

Frank

 

[bigfooted]

what method do you use to actually solve your system? And why do you solve the equations iteratively in a segregated fashion? You could just solve the coupled finite difference equations using a tridiagonal solver (TDMA).

 

[FrankST]

The constant coefficients are sometimes in a way that makes the system singular. So in most cases if I want to solve the system with no iteration (in the sense I mentioned in the first post) I end up with a singular system of equations.

Let me rewrite the system in in actual forms:

1- a1 * f,xx + a2 * f,yy + a3 * g,xx + a4 * g,yy - a5 * f - a6 * g = 0

2- a1 * f,xx + a2 * f,yy + a3 * g,xx + a4 * g,yy - a5 * f - a6 * g = 0

As you see now the constants are the same in the first and the second PDE.

The only thing that makes a distinction between f(x,y) and g(x,y) is their boundary conditions that are different within the domain.

 

[meldraft]

Can you tell us the physical problem that this equation describes? It looks like a linear combination of non-homogeneous Laplace's equations but I can't tell what the phenomenon is.

 

[FrankST]

It is related to the deformation that two indentors make in an elastic membrane.

 

[Chestermiller]

What is the specific problem, and the elastic analysis of the problem?

 

[FrankST]

I've already explained the problem. Do you have any idea on another method to solve the equations?

 

[meldraft]

How about this. Let's say you use a central fd scheme. Now:

$$\frac{\partial^2 f}{\partial x^2}=\frac{f(x+h,y)-2f(x,y)+f(x-h,y)}{h^2}$$

$$\frac{\partial^2 f}{\partial y^2}=\frac{f(x,y+h)-2f(x,y)+f(x,y-h)}{h^2}$$

$$\frac{\partial^2 g}{\partial x^2}=\frac{g(x+h,y)-2g(x,y)+g(x-h,y)}{h^2}$$

$$\frac{\partial^2 g}{\partial y^2}=\frac{g(x,y+h)-2g(x,y)+g(x,y-h)}{h^2}$$

An equivalent problem to yours would be

$$(a_1f,xx+a_2f,yy+a_3g,xx+a_4g,yy-a_5f-a_6g)^2+(b_1f,xx+b_2f,yy+b_3g,xx+b_4g,yy-b_5f-b_6g)^2=0$$

replacing the fd formulas you have an expression where you need to calculate f,g a total of 10 times.

Either way, no matter how you express the problem, the way you solve it of course is starting at the boundaries using a forward fd and going inwards using the central fd. This will usually produce a tridiagonal system, which can be solved in a number of ways. If in your case the resulting system is not tridiagonal, try using a fd scheme that does produce a tridiagonal system. This is how we solve the heat equation. Now what that scheme would be I cannot tell, you will need to try it out for yourself, unless there is someone who has actually solved this before.

 

[FrankST]

First of all, thanks a lot for your attention to this problem.

Back to the problem,

You squared each equation and you summed them up.

Then you used FD equivalent to discretize the equations.

What I can't understand is that after replacing FD formulas inside the expression with "squared" terms you must end up with a nonlinear expression. Am I wrong? (maybe I'm missing a very obvious thing)

If you help me to understand this part then I can comment on the tridiagonal matrix that you mentioned.

 

[FrankST]

Maybe this is what you're meaning:

$(a_1f,xx+a_2f,yy+a_3g,xx+a_4g,yy-a_5f-a_6g)^2+(b_1f,xx+b_2f,yy+b_3g,xx+b_4g,yy-b_5f-b_6g)^2=0$

$\Rightarrow$

$(a_1f,xx+a_2f,yy+a_3g,xx+a_4g,yy-a_5f-a_6g)=+(b_1f,xx+b_2f,yy+b_3g,xx+b_4g,yy-b_5f-b_6g)$

and

$(a_1f,xx+a_2f,yy+a_3g,xx+a_4g,yy-a_5f-a_6g)=-(b_1f,xx+b_2f,yy+b_3g,xx+b_4g,yy-b_5f-b_6g)$

am I right?

 

[meldraft]

I was thinking that you could use the squared expression in order to use optimization methods to achieve faster convergence, since depending on the problem it can converge quite fast. However, before trying that, there might be a couple more things you could try for better convergence.

The first is called relaxation, which consists of the following steps:

1. Set each node inside the boundary to 0 or a convenient value
2. Loop systematically through the interior mesh points, setting each interior point to the average of its four neighbors.
3. Continue this process until the solution converges to the desired accuracy.

This is used for the Laplace equation, but it should be applicable here.

The second thing I had in mind is how you order your nodes. The ordering of the mesh may be responsible for giving you a singular matrix. Try building the equations in a different order.

 

[FrankST]

About the relaxation method:

I can't get the step 2. After assuming some value for the nodes in step 1 I would require to solve Eq1 for f1. Then when should I take the average of the surrounding nodes?

What I can think is that after solving Eq1 for f1 I should update the nodal values for f1 by taking the average of the neighboring nodes and after updating f1 I can substitute f1 in Eq2 and solve Eq2 for f2 and iterate until convergence. Is it what you are suggesting?

I will try your second comment on node ordering and let you the result. However, since I am dealing with thousands of nodes if I want to simultaneously solve for f1 and f2, instead of dealing with NxN matrix equation I will deal with 2Nx2N system which is a bit more expensive, but it is a secondary issue here.

 

[meldraft]

I'm familiar with the Laplace equation implementation, there should be some calculations for this to work here. Check these out:

http://en.wikipedia.org/wiki/Relaxation_(iterative_method)#Model_problem_of_potential_theory

http://www.damtp.cam.ac.uk/lab/people/sd/lectures/nummeth98/pdes.htm (go to 'mathematics of relaxation' and 'Multigrid')

http://www.math.ohiou.edu/courses/math344/lecture34.pdf

You would have to create the proper relaxation scheme for yourself to use the method. Note that this method is not particularly fast, but you can use it as a benchmark to see how it compares to your own. For better convergence you can use the multigrid method which is also described in the second link. I'm afraid though that I may be taking you a little off topic with these methods, so feel free to say so if that is the case. If you want to try them out, you can find online pretty decent material on how to set them up.