4
$\begingroup$

I'm trying to solve with Mathematica the following problem $$-\Delta u = 10$$ on $[0,1]\times [0,1]$ with homogeneous Dirichlet boundary conditions.

I obtain the following plot with my own Matlab code, but with Mathematica I can't even obtain the plot. Could anyone confirm if the same plot is obtained with your own version of Mathematica?

f = 10;
pde = -Laplacian[u[x, y], {x, y}] == f;
bc = {u[x, 0]  == 0,
   u[x,1]  == 0,
   u[0, y] == 0,
   u[1, y] == 0
   };
DSolve[{pde, bc}, u[x, y], {x, y}, Assumptions -> {0 <= x <= 1 && 0 <= y <= 1}]

enter image description here

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Michael E2
    Commented Mar 9, 2021 at 3:21

2 Answers 2

Reset to default
4
$\begingroup$

I'm sure you were doing this numerically in Matlab, so here's the equivalent in Mathematica

whee = NDSolveValue[{pde, bc}, u[x, y],
  {x, 0, 1},
  {y, 0, 1}
  ];
Plot3D[whee, {x, 0, 1}, {y, 0, 1}]

enter image description here

To obtain the max

NMaximize[whee,  {x, y} \[Element] Rectangle[]]

{0.736711, {x -> 0.500081, y -> 0.499919}}
Improve this answer
$\endgroup$
8
  • $\begingroup$ Did you set $f=10$ in your code? @b3m2a1 $\endgroup$
    – bobinthebox
    Commented Mar 9, 2021 at 0:19
  • $\begingroup$ @bobinthebox I just copy-pasted what was there $\endgroup$
    – b3m2a1
    Commented Mar 9, 2021 at 0:20
  • $\begingroup$ So, could you confirm that the commands in my Mathematica snippet are the right way to solve this simple problem? @b3m2a1 ? $\endgroup$
    – bobinthebox
    Commented Mar 9, 2021 at 0:21
  • $\begingroup$ @bobinthebox except use NDSolve and NDSolveValue so you get a numerical solution... $\endgroup$
    – b3m2a1
    Commented Mar 9, 2021 at 0:22
  • $\begingroup$ Thanks for your time, I see now the problem. Before accepting your answer, I have a last curiosity: how can you see the maximum of the solution? For instance, with my own Matlab code, I obtain a maximum of $0.73$. @b3m2a1 $\endgroup$
    – bobinthebox
    Commented Mar 9, 2021 at 0:24
4
$\begingroup$

For some reason, I have to multiply the pde by $-1$ in order to get any results, Mathematica is just idiosyncratic. Consider

f=10;
pde=Laplacian[u[x,y],{x,y}]==-f;
bc={u[x,0]==0,u[x,1]==0,u[0,y]==0,u[1,y]==0};
soln=DSolve[{pde,bc},u[x,y],{x,y}]

The result is a list of a list of rules, containing a rule which defines $u(x,y)$. You can use StringReplace[ToString[TeXForm@u[x,y]/.First@soln],{"K[1]"->"k"}] to obtain $$2 \underset{k=1}{\overset{\infty }{\sum }}-\frac{40 \text{sech}\left(\frac{1}{2} \pi k\right) \sin ^2\left(\frac{1}{2} \pi k\right) \sin(\pi x k) \sinh\left(\frac{1}{2} \pi (y-1) k\right) \sinh\left(\frac{1}{2} \pi y k\right)}{\pi ^3 k^3}$$ We can make an approximate plot with

Plot3D[Sum[-80/(\[Pi]^3 k^3)
Sech[(\[Pi] k)/2] Sin[(\[Pi] k)/2]^2
Sin[\[Pi] x k] Sinh[(\[Pi] (y - 1) k)/2] Sinh[(\[Pi] y k)/2],
{k, 1, 20}], {x, 0, 1}, {y, 0, 1}, PlotStyle -> Green]

soln

Improve this answer
$\endgroup$
5
  • $\begingroup$ It seems pretty close to mine. Could you tell me what's the maximum of the solution ? @Adam $\endgroup$
    – bobinthebox
    Commented Mar 9, 2021 at 0:22
  • 1
    $\begingroup$ My guess is the idiosyncrasy comes from the fact that DSolve has to reduce its PDE argument down to something that can be structurally compared to a set of templates it has built in. I'm sure the -Laplacian is just messing with that $\endgroup$
    – b3m2a1
    Commented Mar 9, 2021 at 0:31
  • $\begingroup$ It's not too hard to show $\partial_xu(1/2,1/2)=\partial_yu(1/2,1/2)=0$ with Mathematica, which shows $(1/2,1/2)$ is a local max. Showing its a global max is harder; not sure if Mathematica can do this without manual work. $\endgroup$
    – Adam
    Commented Mar 9, 2021 at 0:44
  • $\begingroup$ @b3m2a1 Laplacian seems to work ok; even DSolving with -D[D[u[x,y],x],x]-D[D[u[x,y],y],y]==f fails (while D[D[u[x,y],x],x]+D[D[u[x,y],y],y]==-f works just as well as Laplacian) $\endgroup$
    – Adam
    Commented Mar 9, 2021 at 0:46
  • 1
    $\begingroup$ @Adam my comment isn't about Laplacian, per se, but about putting the negative of the differential operator on the left in this specific case. As I said, DSolve compares structurally against a set of templates and it clearly isn't ready for that one. I'd report it to WRI $\endgroup$
    – b3m2a1
    Commented Mar 9, 2021 at 0:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.