Drexel dragonThe Math ForumDonate to the Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by Drexel University or The Math Forum.


Math Forum » Discussions » Software » comp.soft-sys.matlab

Topic: pdenonlin "Stepsize too small"
Replies: 2   Last Post: Jul 9, 2013 10:18 AM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Bill Greene

Posts: 26
Registered: 2/27/12
Re: pdenonlin "Stepsize too small"
Posted: Jul 8, 2013 3:22 PM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply

"Steven Finch" wrote in message <krejde$hde$1@newscl01ah.mathworks.com>...
> Please download two files ellipseb2.m & ellipseg.m from
>
> http://www.people.fas.harvard.edu/~sfinch/matlab/
>
> and execute the following code:
>
> g='ellipseg';
> b='ellipseb2';
> c='1./sqrt(1+ux.^2+uy.^2)';
> a=0;
> f=0;
> rtol=1e-4;
>
> [p,e,t]=initmesh(g);
> [p,e,t]=refinemesh(g,p,e,t);
>
> u=pdenonlin(b,p,e,t,c,a,f,'tol',rtol);
> figure;
> pdeplot(p,e,t,'xydata',u,'zdata',u,'mesh','on','colorbar','off');
> axis equal tight;
>
> format long
>
> [ux,uy]=pdegrad(p,t,u);
> tArea = pdetrg(p,t);
> surfaceArea = sqrt(1 + ux.^2 + uy.^2)*tArea';
> surfaceArea
>
> [p,e,t]=refinemesh(g,p,e,t);
> u=pdenonlin(b,p,e,t,c,a,f,'tol',rtol);
> [ux,uy]=pdegrad(p,t,u);
> tArea = pdetrg(p,t);
> surfaceArea = sqrt(1 + ux.^2 + uy.^2)*tArea';
> surfaceArea
>
> Everything works fine up to this point. The graph shows the minimal surface spanning a folded circular loop (with sharp corners at two boundary points). A good numerical estimate 2.48201... is found (the true value is 2.48228...). If, however, we attempt to improve upon the estimate via
>
> [p,e,t]=refinemesh(g,p,e,t);
> u=pdenonlin(b,p,e,t,c,a,f,'tol',rtol);
>
> then the following message appears:
>
> ??? Error using ==> pdenonlin at 251
> Stepsize too small.
>
> Does anyone have suggestions for how to avoid this problem? Thank you,
>
> Steve Finch
>
> P.S. The two files were created using the command:
>
> pdeellip(0,0,1,1/sqrt(2),0)
>
> then specifying u = abs(y) using the Boundary Mode menu, and finally using the commands:
>
> fid = wgeom(g, 'ellipseg');
> fid = wbound(b, 'ellipseb2');


Steven,

If you add the option:
'Jacobian', 'full'

to your pdenonlin calls, you will get converged solutions with only a few
iterations. I also recommend adding the option, 'Report', 'on'
so you can see what is happening during the Newton-Raphson algorithm.

u=pdenonlin(b,p,e,t,c,a,f,'tol',rtol, 'Jacobian', 'full', 'Report', 'on');

I tried your example in both MATLAB R2012b and R2013a. Convergence is
slightly faster in R2013a due to some improvements made to pdenonlin
in that release.

Bill



Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© Drexel University 1994-2014. All Rights Reserved.
The Math Forum is a research and educational enterprise of the Drexel University School of Education.