finite difference method for nonlinear ode

Question

I am trying to use finite difference method to solve $y'''+y^2y''-y'=0, y(0)=y'(0)=0, y''(1)=1$. I let $u=dy/dx$ so the new problem is $u''+y^2u'-u=0, u(0)=0, u'(1)=1,y=\int u$. To try and solve this in MATLAB I first solved a similar linear problem $u''+u'-u=0$. Following code sets up the two matrices so that we can solve $A=b$

b=zeros(n,1);
b(1)=0;

A=zeros(n,n);
A(1,1)=2+h^2;
A(1,2)=-1-(h/2);

for i=2:n-1
A(i,i-1)=-1+(h/2);
A(i,i)=2+h^2;
A(i,i+1)=-1-(h/2);
end

b(n)=(1+(h/2))*beta;
A(n,n)=2+h^2;
A(n,n-1)=-1+(h/2);

I am having trouble adding the $y$ term to the $A$ matrix. How can I do this? I am planning on using the Simpsons method to find the $y^2$ term.

Answer 1

First, you should split $$y'''+y^2y''-y'=0$$ into a system of three first order ordinary differential equations by letting $$y'= u$$ and $$ u'=v $$ such that you get $$v'+y^2v-u = 0$$

After taking another look, I think I might now see what you are trying to do with Simpson's rule. Starting with the definition of Simpson's rule $$\int_{x}^{x+2\Delta x}f(x)dx=\frac{\Delta x}{3}(f(x)+4f(x+\Delta x)+f(x+2\Delta x))$$ if you consider a grid of $n$ points evenly spaced by $\Delta x$ where $x_i=i \Delta x$ for $i=0,1,...,n-1$, then for your ODE system you have $$y_{i+2}-y_i = \frac{\Delta x}{3}(u_{i}+4u_{i+1}+u_{i+2})$$ $$u_{i+2}-u_i = \frac{\Delta x}{3}(v_{i}+4v_{i+1}+v_{i+2})$$ $$v_{i+2}-v_i = \frac{\Delta x}{3}(u_{i}-y_{i}^2v_{i}+4(u_{i+1}-y_{i+1}^2v_{i+1})+u_{i+2}-y_{i+2}^2v_{i+2})$$

This is not a linear system of equations in $y$, $u$, and $v$. You can still proceed, but it is difficult. Alternatively, you can integrate these equation in any number of other ways. I recommend Matlabs ODE solvers such as Ode45. Maybe I still don't understand what you are referring to with Simpson's method. If you have simple initial conditions such as y(0)=u(0)=v(0)=1 it would be straight forward to use Ode45. However, it appears you have boundary conditions y(0)=u(0)=v(1)=1 which means you will need to use something like a shooting method (root finding) to figure out which inital conditions y(0)=u(0)=v(0)=α lead you to v(1)=1.

Answer 2

The usual way to treat non-linear ODE is to use iterations to resolve the non-linearity.

Starting from the change of variable $u = y'$, you have indeed

$$u'' + y^2u' -u = 0$$

with the conditions $u(0)=0$ and $u'(1)=1$. To start the iterative method, take a (reasonable) initial value for $y$, $y_{0}$. Then, solve for $u_{1}$ the ODE $$u_{1}'' + y_0^2u_1' -u_1 = 0$$

From $u_1$, you can recover $y_1$ using $y_1(x) = \int_{0}^x u_1(s) ds$ (since y(0)=0).

Now, you can go computing $u_2$ from $y_1$, the $y_2$ from $u_2$ and so on, until you're happy with the solution.