## Coupled PDE System - Numerical Solution

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.

Frank
 PhysOrg.com science news on PhysOrg.com >> Galaxies fed by funnels of fuel>> The better to see you with: Scientists build record-setting metamaterial flat lens>> Google eyes emerging markets networks
 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
 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.

## Coupled PDE System - Numerical Solution

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.

 Quote by FrankST 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

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.
 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.
 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.
 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
 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).
 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.
 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.
 It is related to the deformation that two indentors make in an elastic membrane.
 Recognitions: Gold Member What is the specific problem, and the elastic analysis of the problem?
 I've already explained the problem. Do you have any idea on another method to solve the equations?
 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.
 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.
 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?
 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.

 Tags coupled pde, finite difference