How would I solve y'' = sinh(y) using finite differences?

Hello,

I've tried a few things like:

N = 50;

A = toeplitz([2 -1 zeros(1, N-2)]);

x = linspace(-1, 1, N)';

y = A\sinh(x);

Of course, what happens is that is solves for y'' = sinh(x), not sinh(y). How can I make this solve for sinh(y)?

Thanks,

Mark

[Matt Tearle]

Is there a reason it has to be done with FD? Is this a problem for coursework? In which case, are there any further guidelines (regarding what you can/can't/must/should do)?

[Matt Fig]

In addition to Matt's comments, what are the boundary conditions?

[Mark]

Matt, this isn't a problem for coursework. I don't have to use finite differences, but I would like to solve it numerically as a boundary value problem.

The boundaries I'm interested in are generally y(C1) = C2, y(C3) = C4 where C1, C2, C3, and C4 are constants.

Thanks.

[Mark]

I'm doing this as part of my PhD work, unfortunately differential equations (and math in general) is not my strong point.

Answer 1

[Matt Tearle]

Great, if it doesn't have to be a specific FD formulation, you can use bvp4c. That link shows an example of a second-order equation.

But if you really want to do it yourself, you should have access to Optimization Toolbox, so you could use fsolve to solve the system of equations F(y) = A*y - sinh(y) = 0 for the vector y.

Or, if you really really want to go it alone, you could implement a scheme to solve that system, such as a Newton method. The Jacobian is just A - diag(cosh(y)). Then iterate y_[n+1] = y_n - dy where dy = J\F.

EDIT: and because I am a massive dork:

c1 = 2;  c2 = pi;  c3 = 7;  c4 = 4.2;

xi = linspace(c1,c3);
yi = linspace(c2,c4);
dy = @(t,y) [y(2);sinh(y(1))];
bc = @(ya,yb) [ya(1)-c2;yb(1)-c4];

ysol = bvp4c(dy,bc,bvpinit(xi,[0,1]));

y = deval(ysol,xi);
plot(xi,y(1,:))

More EDITs I cleaned up a typo, but also... note that the discretization stepsize is not accounted for. A should actually be replaced by A/h^2 and diag(cosh(y)) by diag(cosh(y)/h). Also, the boundary conditions are not accounted for, either. Finally, the A Mark defined is the negative of the usual definition of y'', so either flip the sign there or in the definition of F and J.

[Mark]

Rock on! Thanks a lot that's just what I was looking for.