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:

[tex]N(x,y)=ex1^2+ex2^2[/tex]

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.

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

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:

[tex]\frac{\partial^2 f}{\partial x^2}=\frac{f(x+h,y)-2f(x,y)+f(x-h,y)}{h^2}[/tex]

[tex]\frac{\partial^2 f}{\partial y^2}=\frac{f(x,y+h)-2f(x,y)+f(x,y-h)}{h^2}[/tex]

[tex]\frac{\partial^2 g}{\partial x^2}=\frac{g(x+h,y)-2g(x,y)+g(x-h,y)}{h^2}[/tex]

[tex]\frac{\partial^2 g}{\partial y^2}=\frac{g(x,y+h)-2g(x,y)+g(x,y-h)}{h^2}[/tex]

An equivalent problem to yours would be

[tex](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[/tex]

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:

[itex](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[/itex]

[itex]\Rightarrow[/itex]

[itex](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)[/itex]

and

[itex](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)[/itex]

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.