Using MATLAB to solve Poisson matrix equation

2.6k Views Asked by At

I'm not sure if this is the correct forum to post my question. I'm working on a Poisson-based maths assignment and am stuck as regards finding the solution to the Poisson matrix equation. The matrix I have is rectangular. I've found some MATLAB code online for solving Poisson's equation and am just wondering if you could suggest which might be the best to look into for my particular problem (question 4)?

Poisson questions

I'm new to Poisson and MATLAB, so thanks for any advice you can give.

Links to MATLAB code for solving Poisson's equation:

https://uk.mathworks.com/matlabcentral/fileexchange/38090-2d-poisson-equation

https://uk.mathworks.com/help/pde/ug/pde.pdemodel.solvepde.html

https://staff.washington.edu/rjl/fdmbook/matlab/poisson.m

https://www.math.utah.edu/~wright/courses/5620/notes/FD2PoissonEx.pdf

http://people.bu.edu/andasari/courses/Fall2015/LectureNotes/Lecture14_27Oct2015.pdf

https://math.boisestate.edu/~wright/courses/m567/FD2-Poisson.pdf

Answers 1 and 2

Answer 3 Answer 3 Answer 3 Answer 3 Answer 3

MATLAB output:

ans =

[ -14203947/3248336,  -940473/203021, -382635/141232,     233/203021,  6708477/3248336]
[     -362575/31234,   -217820/15617,    -16915/1358,  -146040/15617,    -171485/31234]
[ -18233133/3248336, -1440217/203021, -891885/141232, -905553/203021, -7065717/3248336]

$\therefore$

$\begin{align} P_{13} = -4.37\qquad P_{23} = -4.63\qquad P_{33} = -2.71\qquad P_{43} = 0.00\qquad\,\,\,\,\, P_{53} = 2.07 \end{align}$

$\begin{align} P_{12} = -11.61\,\,\,\,\,\,\,\,\, P_{22} = -13.95\,\,\,\,\,\,\,\,\, P_{32} = -12.46\,\,\,\,\,\,\,\,\, P_{42} = -9.35\qquad P_{52} = -5.49 \end{align}$

$\begin{align} P_{11} = -5.61\qquad P_{21} = -7.09\qquad P_{31} = -6.32\qquad P_{41} = -4.46\qquad P_{51} = -2.18 \end{align}$

1

There are 1 best solutions below

2
On BEST ANSWER

Since you have shown your work, I decide to write an answer. (Although I think your question is more on-topic for stackoverflow. Also, please write in latex instead of posting images in the future since images are impossible to search (our OCR technology is not that good yet))

Actually, you have already solved the harder part of the problem. Assuming you've done it right, the only thing left to do is to solve the linear system given by the system of equations: $$P_{i, j} = \frac{P_{i + 1, j} + P_{i - 1, j} + P_{i, j + 1} + P_{i, j - 1} - 25 + 5 \frac{i}{2} \frac{j}{2}}{4}$$ where $1\le i \le 5$ and $1 \le j \le 3$ with the convention that $P_{i, j} = \begin{cases} 20 &\text{ if } j < 1 \text{ or } j > 3 \\ 0 &\text{ if (} i < 1 \text{ or } i > 5 \text{) and } 1 \le j \le n \end{cases}$ (Thank god we don't have buffer overflow!)

The following program does exactly that. I have tested it in octave with octsympy. Please check if it works in matlab as well.

%{ see https://stackoverflow.com/questions/35115039 %}
if__ = @(varargin) varargin{3 - varargin{1}}();
if_ = @(pred, if_clause, else_clause) if__(pred, if_clause, else_clause);

%{ size of the matrix %}
m = 5;
n = 3;

%{ if j goes out of range, return 20 %}
%{ else if i goes out of range, return 0 %}
%{ otherwise, return the symbol Pij %}
P = @(i, j) if_(j < 1 || j > n, ...
                @() 20, ...
                @() if_(i < 1 || i > m, ...
                        @() 0, ...
                        @() sym(sprintf('P%d%d', i, j))));

%{ our equation with P(i, j) moved to the RHS %}
eqn = @(i, j) ...
       (P(i + 1, j) + P(i - 1, j) + P(i, j + 1) + P(i, j - 1) ...
        - 25 ...
        + sym(5) * i / 2 * j / 2) ...
      / 4 ...
      - P(i, j);

%{ system of equations obtained by substituting the above equation %}
%{ with 1 <= i <= m and 1 <= j <= n %}
eqns = arrayfun(@(i) arrayfun(@(j) eqn(i, j), ...
                              1:n, ...
                              'UniformOutput', false), ...
                1:m, ...
                'UniformOutput', false);
eqns2 = vertcat(eqns{:});
eqns3 = vertcat(eqns2{:});

%{ solve the system of equations %}
solns = solve(eqns3);

%{ display the solutions %}
[solns.P13 solns.P23 solns.P33 solns.P43 solns.P53; ...
 solns.P12 solns.P22 solns.P32 solns.P42 solns.P52; ...
 solns.P11 solns.P21 solns.P31 solns.P41 solns.P51]