1
$\begingroup$

I am trying to solve the differential equation $y''(x) = y(x)+\sin[y'(x)]$ numerically over the interval $[0,1].$ Now using central difference method I arrived at the equations of the type $$ y_k = \frac{1}{h^2+2} \left(y_{k-1} + y_{k-1} - h^2\sin\left(\frac{y_{k+1}-y_{k-1}}{2h}\right)\right) $$

for $k=1,2,\ldots, m-1.$ I am trying to solve for $m=10$ making a system of equations with $y_0=0$ and $y_{10}=1.$ The issue is I am not able to produce the equations in Mathematica.

SysEq[xs_List] := (1/(h^2 + 2) (#[[1]] + #[[3]] - h^2*Sin[(#[[3]] - #[[1]])/(2*h)])) & /@ Partition[xs, 10, 1, {1, 1}]; 

So that I can call it

SysEq[{x1, x2, x3, x4, x5, x6, x7, x8, x9, x10}]

or may be an easy way to create whole system of equations.

$\endgroup$

1 Answer 1

2
$\begingroup$

to generate the equations, just use Table command. To solve them is different issue. Since they are non-linear.

$$ y^{\prime\prime}=y+\sin\left( y^{\prime}\right) $$

3 point centered difference for $y^{\prime\prime}$ is

$$ y^{\prime\prime}=\frac{1}{h^{2}}\left( y_{n-1}-2y_{n}+y_{n+1}\right) $$

centered difference for $y^{\prime}$ is

$$ y^{\prime}=\frac{1}{2h}\left( y_{n+1}-y_{n-1}\right) $$

The ode becomes

$$ \frac{1}{h^{2}}\left( y_{n-1}-2y_{n}+y_{n+1}\right) =y_{n}+\sin\left( \frac{1}{2h}\left( y_{n+1}-y_{n-1}\right) \right) $$

let $L=1$ and 11 points. 9 are internal points, and 2 are boundary points. Hence number of intervals is 10.

Therefore $h=\frac{1}{10}$ and $n=0,1,\cdots,10$. Where $y_{0}=0,y_{10}=1$. For $n=1$

$$ \frac{1}{h^{2}}\left( y_{0}-2y_{1}+y_{2}\right) =y_{1}+\sin\left( \frac {1}{2h}\left( y_{2}-y_{0}\right) \right) $$

For $n=2$

$$ \frac{1}{h^{2}}\left( y_{1}-2y_{2}+y_{3}\right) =y_{2}+\sin\left( \frac {1}{2h}\left( y_{3}-y_{1}\right) \right) $$

And so on. For $n=9$

$$ \frac{1}{h^{2}}\left( y_{8}-2y_{9}+y_{10}\right) =y_{9}+\sin\left( \frac {1}{2h}\left( y_{10}-y_{8}\right) \right) $$

Hence the equations are

\begin{align*} y_{0} & =0\\ \frac{1}{h^{2}}\left( y_{0}-2y_{1}+y_{2}\right) & =y_{1}+\sin\left( \frac{1}{2h}\left( y_{2}-y_{0}\right) \right) \\ \frac{1}{h^{2}}\left( y_{1}-2y_{2}+y_{3}\right) & =y_{2}+\sin\left( \frac{1}{2h}\left( y_{3}-y_{1}\right) \right) \\ & \vdots\\ \frac{1}{h^{2}}\left( y_{8}-2y_{9}+y_{10}\right) & =y_{9}+\sin\left( \frac{1}{2h}\left( y_{10}-y_{8}\right) \right) \\ y_{10} & =1 \end{align*}

Code is

numberInternalPoints = 9;
h                    = 1/(numberInternalPoints+1); 
y[0]                 = 0; 
y[numberInternalPoints + 1] = 1;

vars = Table[y[n], {n, 1, numberInternalPoints}]
eqs = Table[
   1/h^2*(y[n - 1] - 2 y[n] + y[n + 1]) == 
    y[n] + Sin[1/(2*h)*(y[n + 1] - y[n - 1])], {n,1, numberInternalPoints}];
TableForm[eqs]

enter image description here

9 equations, 9 unknowns.

Update

As suggested by user64494 below, FindRoot can be used to solve these equations, and the result agree with NDSolve

sol = FindRoot[
  eqs, {{y[1], 0}, {y[2], 0}, {y[3], 0}, {y[4], 0}, {y[5], 0}, {y[6], 
    0}, {y[7], 0}, {y[8], 0}, {y[9], 0}}]
data = Table[{n*h, Last@sol[[n]]}, {n, 1, numberInternalPoints}]
data = PrependTo[data, {0, 0}];
data = AppendTo[data, {1, 1}]

Gives

{{0, 0}, {1/10, 0.0558825}, {1/5, 0.117883}, {3/10, 0.187164}, {2/5, 
  0.265029}, {1/2, 0.352914}, {3/5, 0.452386}, {7/10, 0.565109}, {4/5,
   0.692812}, {9/10, 0.837223}, {1, 1}}

And

p = ListLinePlot[data, Mesh -> All, MeshStyle -> Red]
ClearAll[y, x]
sol = NDSolveValue[{y''[x] == y[x] + Sin[y'[x]], y[0] == 0, 
    y[1] == 1}, y, {x, 0, 1}];
p2 = Plot[sol[x], {x, 0, 1}, PlotStyle -> Blue];
Show[p, p2]

enter image description here

$\endgroup$
6
  • $\begingroup$ A spot on the Sun: NSolve[eqs, vars, Reals] generates "NSolve::nsmet: This system cannot be solved with the methods available to NSolve". $\endgroup$
    – user64494
    Commented Sep 21, 2023 at 19:21
  • $\begingroup$ Just to compare. The command of Maple DirectSearch:-SolveEquations([100*(-2*y[1]+y[2]) = sin(5*y[2])+y[1], 100*(y[1]-2*y[2]+y[3]) = sin(5*(-y[1]+y[3]))+y[2], 100*(y[2]-2*y[3]+y[4]) = sin(5*(-y[2]+y[4]))+y[3], 100*(y[3]-2*y[4]+y[5]) = sin(5*(-y[3]+y[5]))+y[4], 100*(y[4]-2*y[5]+y[6]) = sin(5*(-y[4]+y[6]))+y[5], 100*(y[5]-2*y[6]+y[7]) = sin(5*(-y[5]+y[7]))+y[6], 100*(y[6]-2*y[7]+y[8]) = sin(5*(-y[6]+y[8]))+y[7], 100*(y[7]-2*y[8]+y[9]) = sin(5*(-y[7]+y[9]))+y[8], 100*(1+y[8]-2*y[9]) = sin(5*(1-y[8]))+y[9]]) results in $\endgroup$
    – user64494
    Commented Sep 21, 2023 at 19:46
  • $\begingroup$ [y[1] = 0.558824815612752e-1, y[2] = .117882524237100, y[3] = .187164143110649, y[4] = .265028696216958, y[5] = .352914417317458, y[6] = .452385864615537, y[7] = .565109488648346, y[8] = .692812302307770, y[9] = .837223071943341] in a moment. $\endgroup$
    – user64494
    Commented Sep 21, 2023 at 19:47
  • 1
    $\begingroup$ FindRoot[{100 (-2 y[1] + y[2]) == Sin[5 y[2]] + y[1], 100 (y[1] - 2 y[2] + y[3]) == Sin[5 (-y[1] + y[3])] + y[2], 100 (y[2] - 2 y[3] + y[4]) == Sin[5 (-y[2] + y[4])] + y[3], 100 (y[3] - 2 y[4] + y[5]) == Sin[5 (-y[3] + y[5])] + y[4], 100 (y[4] - 2 y[5] + y[6]) == Sin[5 (-y[4] + y[6])] + y[5], 100 (y[5] - 2 y[6] + y[7]) == Sin[5 (-y[5] + y[7])] + y[6], 100 (y[6] - 2 y[7] + y[8]) == Sin[5 (-y[6] + y[8])] + y[7], 100 (y[7] - 2 y[8] + y[9]) == Sin[5 (-y[7] + y[9])] + y[8], 100 (1 + y[8] - 2 y[9]) == Sin[5 (1 - y[8])] + y[9]}, $\endgroup$
    – user64494
    Commented Sep 21, 2023 at 20:03
  • $\begingroup$ {{y[1], 0}, {y[2], 0}, {y[3], 0}, {y[4], 0}, {y[5], 0}, {y[6], 0}, {y[7], 0}, {y[8], 0}, {y[9], 0}}] outputs {y[1] -> 0.0558825, y[2] -> 0.117883, y[3] -> 0.187164, y[4] -> 0.265029, y[5] -> 0.352914, y[6] -> 0.452386, y[7] -> 0.565109, y[8] -> 0.692812, y[9] -> 0.837223} $\endgroup$
    – user64494
    Commented Sep 21, 2023 at 20:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.