Sorry for the confusion, the numerics for linear ode solving in ApproxFun is built on the ultra spherical spectral method by myself and Alex Townsend:
A fast and well-conditioned spectral method, SIAM Review, 55: 462–489 http://epubs.siam.org/doi/abs/10.1137/120865458 It's different from existing bvp approaches for several reasons: it is stable and O(n) with guaranteed accuracy. The efficient O(n) implementation of the approach is one of the reasons I got interested in Julia: I implemented originally in C++ but the data structure manipulation was not conducive to research Sent from my iPad > On 21 Jun 2014, at 10:58 pm, 'Stéphane Laurent' via julia-users > <[email protected]> wrote: > > Sorry but I really don't understand you first sentence. > About the code, I will possibly experiment it during the next days or weeks, > and I'll send you some feedback. > Thank you again for this library and your help. > - Stéphane > > Le samedi 21 juin 2014 14:24:03 UTC+2, Sheehan Olver a écrit : >> >> Thanks! >> >> I'm interested to know if the ultraspherical spectral approach has any >> clear advantages over spectral collocation for nonlinear odes. (For linear >> odes the advantages are clear.) >> >> So if you find the code better than expected please let me know. It's a bit >> buggy and unoptimized, so chances are you won't. >> >> Sent from my iPad >> >>> On 21 Jun 2014, at 10:15 pm, 'Stéphane Laurent' via julia-users >>> <[email protected]> wrote: >>> >>> Ok, I understand. It works and it is really awesome. Thank you ! >>> >>> >>> Le samedi 21 juin 2014 12:17:22 UTC+2, Sheehan Olver a écrit : >>>> >>>> Hi, >>>> >>>> >>>> You didn’t quite get the Newton iteration right: if you want to solve >>>> >>>> B u + [1,0] = 0 >>>> L u + g(u) = 0 >>>> >>>> then Newton iteration becomes >>>> >>>> u = u - [B, L + g’(u)]\[B u + [1,0], L u + g(u)] >>>> >>>> >>>> i.e., your bc right hand side is not right. Below is the corrected code. >>>> >>>> Cheers, >>>> >>>> Sheehan >>>> >>>> >>>> >>>> x=Fun(identity, Interval(0.,1.)) >>>> d=x.domain >>>> B=neumann(d) >>>> D=diff(d) >>>> # Solves Lu + g(u)-1==0 >>>> L = D^2 + 2/(x.+1)*D >>>> g = u -> -(exp(u)-exp(-u)); gp = u -> -(exp(u)+exp(-u)) >>>> >>>> u=-0.3*x+0.5 #initial guess >>>> >>>> for k=1:5 # this crashes if Ii increase the number of iterations >>>> u=u-[B,L+gp(u)]\[diff(u)[0.]+1.,diff(u)[1.],L*u+g(u)]; >>>> end >>>> >>>> plot(u) >>>> >>>>> On 21 Jun 2014, at 6:54 pm, 'Stéphane Laurent' via julia-users >>>>> <[email protected]> wrote: >>>>> >>>>> Hello, >>>>> >>>>> I'd like to solve this equation with Neumann boundary conditions. My code >>>>> below does not work. Am I doing something bad or is it a failure of the >>>>> Newton algorithm ? >>>>> >>>>> # solves u" = (exp(u)-exp(-u))-2u'/(x+1) , u'(0)=-1, u(1)=0 >>>>> x=Fun(identity, Interval(0.,1.)) >>>>> d=x.domain >>>>> B=neumann(d) >>>>> D=diff(d) >>>>> # Solves Lu + g(u)-1==0 >>>>> L = D^2 + 2/(x.+1)*D >>>>> g = u -> -(exp(u)-exp(-u)); gp = u -> -(exp(u)+exp(-u)) >>>>> >>>>> u=-0.3*x+0.5 #initial guess >>>>> >>>>> for k=1:5 # this crashes if Ii increase the number of iterations >>>>> u=u-[B,L+gp(u)]\[-1.,0.,L*u+g(u)]; >>>>> end >>>>> >>>>> >>>>> The solution should look like that : >>>>> >>>>> >>>>> >>>>> >>>>
