Solving differential equations using implicit methods

[Sagekilla]

Hi all, I'm writing myself a ordinary differential equation solver and I've already implemented several explicit integrators, which were pretty easy for me to do. Now I've decided to work on some implicit methods (for any stiff equations) and I've run into some issues.

The most basic one is the Backward Euler:
$$y_{n+1} = y_{n} + hf(t_{n+1}, y_{n+1})$$

And as I understand, to solve this you need to use some iteration to get the value of $$y_{n+1}$$ to solve that equation for the next time step.

The problem I have is: How can I solve this? I thought I could use something like Newton's method, but I don't know how to properly apply it.

Any help, even a push in the right direction would be appreciated. Thanks!

 

[matematikawan]

The backward method
$$y_{n+1} = y_{n} + hf(t_{n+1}, y_{n+1})$$

need to solve for y_n+1 via NR method.

$$y_{n+1} - y_{n} - hf(t_{n+1}, y_{n+1})=0$$

Let $$g(z) = z - y_{n} - hf(t_{n+1}, z)$$

The iteration
$$z_{i+1} = z_i - \frac{g(z_i)}{g'(z_i)}$$
should converges to y_n+1 using the initial estimation z₀=y_n.

 

[Sagekilla]

Thanks for your help! I never understood how to properly apply Newton's method before. I implemented it as you had shown, with:

$$z_{i+1} = z_{i} - \frac{z_{i} - y_{n} - hf(t_{n+1}, z_{i})}{1 - hf'(t_{n+1}, z_{i})}$$

Is it safe to assume that I can say that: $$f'(t_{n+1}, z_{i}) = f(t_{n+1}, f(t_{n+1}, z_{i}))$$

I thought it out for a few minutes, and it seemed reasonable in the context of a few examples I put together:

If we let $$\frac{dy}{dt} = -ky$$ then $$f(t, y) = -ky$$ and $$f'(t, y) = -kf(t, y) = f(t, f(t, y))$$

 

[matematikawan]

I thought that we have this relation: t_n+1=t _n + h
where h is the step size increment.
In the case of $\frac{dy}{dt} = -ky$ , f(t_n+1, z) = -kz

 

[Sagekilla]

Yes, there is. But what I was asking was if I could replace the derivative of f inside of Newton's method by f(t, f(t, y)), since for first order ODEs it would be equivalent to differentiating (at least I think it is) all the y components of differential equation.

Because otherwise if I couldn't do that, then I know I'd have to numerically approximate the derivative as well.

 

[matematikawan]

On a second thought, I think it is easier to use the fixed-point interation
z_n+1=g(z_n)
to approximates y_n+1.

This eliminate the requirement to compute the derivative of f.