I want to define a function $f_c(x,y)$ as the solution of the PDE $\Delta f_c(x,y) = g_c(x,y)$, where $g_c$ is a given function depending on a constant $c$. I tried the following:

domain = Rectangle[{0, 0}, {1, 1}];
solution[c_] := NDSolve[{Laplacian[h[x, y], {x, y}] == x + c*y, DirichletCondition[h[x, y] == 0, True]}, h, {x, y} \[Element] domain];
mem : f[x_?NumericQ, y_?NumericQ, c_?NumericQ] := mem = Part[Evaluate[h[x, y] /. solution[c]], 1];
f[0.1, 0.2, 3]
Plot3D[f[x, y, 3], {x, y} \[Element] domain]

The fourth line returns the approximate value of $f_3(0.1,0.2)$, which to me suggests that everything works. But the plot fails saying NDSolve (...) is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing and a bunch of other similar errors.

I understand the problem probably lies in the way I define $f$, but I don't see how to correct it.


This is what ParametricNDSolve is for

domain = Rectangle[{0, 0}, {1, 1}];

Create a parametric function in the parameter c:

pfun = ParametricNDSolveValue[{Laplacian[h[x, y], {x, y}] == x + c*y, 
    DirichletCondition[h[x, y] == 0, True]}, 
   h, {x, y} \[Element] domain, c];

Use a specific value for c to create the interpolating function:

fun = pfun[3];
Plot3D[fun[x, y], {x, y} \[Element] domain]

enter image description here

  • $\begingroup$ Thanks, but I get problems when I want to take g[x_, y_, c_] := Sum[Part[coefficient, i]*g[x, y, i], {i, 1, c}, with c an integer.. Then I get the error The expression i cannot be used as a part specification when trying to evaluate ParametricNDSolveValue. $\endgroup$ Nov 18 '21 at 8:26
  • $\begingroup$ How is you comment related to the question? $\endgroup$
    – user21
    Nov 19 '21 at 6:22
  • $\begingroup$ The function $g_c$ I want to use, uses the parameter $c$ as a summation limit ($x + cy$ was perhaps a wrong MWE). ParametricNDSolveValue does not appear to be able to handle that, hence my question. Is it worthy of a separate post? $\endgroup$ Nov 19 '21 at 10:51
  • $\begingroup$ @BrazilianCérebro, yes, make new post about that. $\endgroup$
    – user21
    Nov 19 '21 at 14:30

You are making things way too complicated. I could not understand it. I like to keep things simple myself. Why not just do this?

ClearAll[c, x, y, h, sol];

domain = Rectangle[{0, 0}, {1, 1}];

sol[c_?NumericQ] := NDSolve[{Laplacian[h[x, y], {x, y}] == x + c*y, 
    DirichletCondition[h[x, y] == 0, True]}, h, Element[{x, y}, domain]];

Plot3D[Evaluate[h[x, y] /. sol[3]], Element[{x, y}, domain]]

Mathematica graphics


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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