0
$\begingroup$

I am trying to find the first Eigenfunction of the Laplacian (in 1D), i.e. a solution of $$ u''(x)=k u(x)\\ u(0)=u(1)=0 $$ with minimal $k>0$ (in this trivial example, I actually know the analytic solution but this is not the point; neither is using any built-in Eigensolver. I simply want to get acquainted with Mathematica).

To do so, I want to use the Power iteration method:

Repeatedly applying the inverse Laplacian to an arbitrary initial function will produce a sequence converging to the solution.

What I currently have is

ClearAll[h];
ClearAll[f];
f[x_] = 0.5 - Abs[x - 0.5];
steps = 1;
Do[
 s = NDSolve[{h''[x] == f[x], h[0] == 0, h[1] == 0}, h, {x, 0, 1}];
 f[x_] = h[x]/h[0.5] /. s;
 , steps]

However, when I take steps=2 or more, I get

   Dot::rect: Nonrectangular tensor encountered.

Note that this is my very first use of Mathematica, and you might have to explain obvious things to me.

Side question: How do I properly output intermediate results in the loop? If I use

f[x_] = h[x]/h[0.5] /. s; ?f

I do get info on f but if I put

f[x_] = h[x]/h[0.5] /. s; ?f ; ?s

I get

Information::nomatch: "No symbol matching ?s found."

(Adding newlines or using commas doesn't help either. For example,

f[x_] = h[x]/h[0.5] /. s; 
?f ,
?s ,

gives no output at all, whereas

f[x_] = h[x]/h[0.5] /. s; 
 ?f ;
  ?s ;

gives

Information::ssym: "\!\(f; Information[\"s\", LongForm -> False]\) is not a symbol or a valid string pattern"

)

$\endgroup$
1
  • $\begingroup$ So at each step the definition of f and s will be interpolating functions. What information do you want to have about them? How do you want to display it? $\endgroup$
    – Jason B.
    May 11, 2016 at 10:43

2 Answers 2

1
$\begingroup$

To solve the first issue, write

 s = NDSolve[{h''[x] == f[x], h[0] == 0, h[1] == 0}, h, {x, 0, 1}][[1]];

instead (notice the final brackets).

$\endgroup$
1
0
$\begingroup$

This allows you to visualize the (small) changes to the function over the course of the iterations,

ClearAll[h];
ClearAll[f];
f[x_] = 0.5 - Abs[x - 0.5];
steps = 4;
list = Join[
   {{#, f[#]} & /@ Subdivide[1, 100]},
   Table[
    s = NDSolveValue[{h''[x] == f[x], h[0] == 0, h[1] == 0}, 
      h, {x, 0, 1}];
    f[x_] = s[x]/s[0.5];
    {#, f[#]} & /@ Subdivide[1, 100]
    , steps]];
ListLinePlot[list, PlotLegends -> Range[0, 4]]

Mathematica graphics

I've also used NDSolveValue to get around your first problem.

$\endgroup$

Your Answer

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

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