daniel andor

Hi All,

Can anyone help me on how to difference an equation such as:
(y'')^2 + y'(y' + y) + y^3 = 0
where the nonlinearity is in the derivative terms too?

How would I adapt a usual central difference scheme such as
y' ~ (y_+1 - y_-1)/(2h) ?

Do I just substitute in the y' and y'' ******** differences and expand??

Thanks,
Daniel.

arnold neumaier

yes, except that you probably need to use formulas of the same
order for y' and y'', or else rewrite the equation as a first
order system.

But there is ready-made software - see links at
http://www.mat.univie.ac.at/~neum/software.html#dae
[~ is a tilde] - which is probably much more reliable than
what you'd be able to create.

Arnold Neumaier

daniel andor

You mean use 2nd or 4th (etc) order differencing for both? (Assuming for
the moment that I can't rewrite it as a first order system.)


Thanks for that.
D

arnold neumaier

no, to get the same accuracy you need higher order for y''
second order differences for y' are O(h^2) accurate, while
second order differences for y'' are only O(h). You need
at least third order differences for y'' to have an error
of O(h^2).

Arnold Neumaier

spellucci

In article <3FA2A717.1F5F5F01@univie.ac.at>,
Arnold Neumaier <Arnold.Neumaier@univie.ac.at> writes:

?????????????????????????????
we speak about the equidistant case?
(y(x-h)-2y(x)+y(x+h))/h^2 = y''(x) + (1/12)y^(4)(xi) h^2

but ... if we discretize the ode this way, will the discrete system be solvable
at all? this ode has a very bad structure since even the simplest discretization
will not give a reasonable structure in the Jacobian matrix
(nothing M-, diagonally dominant .. )
direct use of higher order fnite differences will make things even worse.
(better use Richardson extrapolation to increase the order)
unfortunately I cannot see how transformation into a first order system
y1=y
y2=y1'
(y2')^2 +y2(y2+y1) + y1^3 =0
would make things much easier.
You said nothing about intial / or boundary values. If you can decide
on the sign of y2', one could write

y1' = y2
y2' = sign * sqrt( -y1^3 + y2*(y2+y1))
for an initial value problem , this could be integrated directly
(with some care considering the argument of sqrt)
and for a boundary value problem collocation using a gauss rule comes into mind.
here even more, due to the nonlinearity and the argment of sqrt, problems
will arise. and for the direct setup by finite differences, you will have
trouble with the Jacobian.
hth
peter

arnold neumaier

Yes, of course, in the equidistant case one has superconvergence; sorry.

But one cannot expect this for an implicit system; this kind of structure is
already lost when discretizing the simple
linear differential equation eps*y''+y'+x=0. There is no way of
getting one of these nice properties here.

Arnold Neumaier

daniel andor

Thanks for the advice.

I'm getting a bit of a headache as to how to construct the matrix to
solve:

If F(y,y',y'',x)=0 is my DAE, then I've constructed a discreet version, as
per the differencing on the equidistant grid you mention before. So now I
have a vector function f(y,x)_i=0 which depends on the solution vector y.
(There is also the small matter of the boundary contitions at f_1 and
f_N.)

I can even construct the Jacobian: J(y,x)_ij = df_i/dy_j, except now J_ij
depends on y as well as x, since f_i is nonlinear in those. This looks
good, and is block diagonal.

I also have Dirichlet boundary conditions at the two ends. How does the
Newton iteration proceed?

Am I missing anything important?

Thanks,
Daniel.

daniel andor

Let me save everyone the trouble and answer myself myself.

I solve J(y,x) . deltay = b - f(y,x)
with b the boundary conditions and y -> y + deltay the iterative step.

It actually works, and converges, as long as I increase the nonlinearity
slowly from zero (where the system is linear) to 1 (where the system is
nonlinear). Otherwise if the process tends to blow up quite badly.

Regards,
Daniel.