Convert 2nd order ODE to system of 1st order

[Haydo]

Homework Statement

Convert the following second-order differential equation into a system of first-order equations and solve y(1) and y'(1) with 4th-order Runge-kutta for h=0.5.
$y''(t)+sin(y(t))=0,\ y(0)=1,\ y'(0)=0$

Homework Equations

The Runge-kutta method might be applicable, but I know how to do that part no problem.

The Attempt at a Solution

The issue I'm having is converting it into a proper system of equations.

I let:
$x_1=y,\ x_2=y'$
Thus:
$x_1'=x_2,\\ x_2'=y''=-sin(y)=-sin(x_1)$
with
$x_1(0)=1,\ x_2(0)=0$

This yields:
$\begin{pmatrix}x_1\\x_2\end{pmatrix}'= \begin{pmatrix}x_2\\-sin(x_1)\end{pmatrix}$

I need to put this in the form:
$\overrightarrow{x'}=A\overrightarrow{x}$
so that I can solve it using runge-kutta. Any ideas?

 

[pasmith]

Homework Statement

Convert the following second-order differential equation into a system of first-order equations and solve y(1) and y'(1) with 4th-order Runge-kutta for h=0.5.
$y''(t)+sin(y(t))=0,\ y(0)=1,\ y'(0)=0$

Homework Equations

The Runge-kutta method might be applicable, but I know how to do that part no problem.

The Attempt at a Solution

The issue I'm having is converting it into a proper system of equations.

I let:
$x_1=y,\ x_2=y'$
Thus:
$x_1'=x_2,\\ x_2'=y''=-sin(y)=-sin(x_1)$
with
$x_1(0)=1,\ x_2(0)=0$

This yields:
$\begin{pmatrix}x_1\\x_2\end{pmatrix}'= \begin{pmatrix}x_2\\-sin(x_1)\end{pmatrix}$

Correct so far.

I need to put this in the form:
$\overrightarrow{x'}=A\overrightarrow{x}$
so that I can solve it using runge-kutta. Any ideas?

The ODE isn't linear with constant coefficients, so it can't necessarily be put into that form. But it isn't necessary to do so to use Runge-Kutta (or any other ODE solution algorithm). As long as you have a first-order vector ODE of the form <br /> \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}) and an initial condition on \mathbf{x} you can proceed.

 

[Haydo]

Ok, so I wrote up the following code,

#System of Equations - Problem 5
#Hayden Smotherman
#Math 428
#30 April 2015

x1=1; #y(0)
x2=0; #y'(0)
x= [x1,x2];
step=0.5; #Step size
boundary = 1;
n=boundary/step;

for i=1:n
  k1=[x2,-sin(x1)];
  k2=[x2+step*k1(1)/2,-sin(x1+step*k1(2)/2)];
  k3=[x2+step*k2(1)/2,-sin(x1+step*k2(2)/2)];
  k4=[x2+step*k3(1),-sin(x1+step*k3(2))];
  x=x+(step/6)*(k1+2*k2+2*k3+k4);
end
x

which gives me x= 1.000000000000000 -0.719893726895089. This feels wrong to me, as x(1)=y(t). Any idea what' s up?

 

[pasmith]

Ok, so I wrote up the following code,
#System of Equations - Problem 5
#Hayden Smotherman
#Math 428
#30 April 2015

x1=1; #y(0)
x2=0; #y'(0)

x= [x1,x2];
step=0.5; #Step size
boundary = 1;
n=boundary/step;

for i=1:n
k1=[x2,-sin(x1)];
k2=[x2+step*k1(1)/2,-sin(x1+step*k1(2)/2)];
k3=[x2+step*k2(1)/2,-sin(x1+step*k2(2)/2)];
k4=[x2+step*k3(1),-sin(x1+step*k3(2))];
x=x+(step/6)*(k1+2*k2+2*k3+k4);
end
x

which gives me x= 1.000000000000000 -0.719893726895089. This feels wrong to me, as x(1)=y(t). Any idea what' s up?

Take a closer look at the assignments to k1, etc. Do you want k1 = [x2,-sin(x1)] or k1 = [x(2),-sin(x(1))]?

 

[Haydo]

Haha, x(2) or it's not going to iterate. Thanks, that was a silly mistake.

 

[D H]

Ok, so I wrote up the following code,

#System of Equations - Problem 5
#Hayden Smotherman
#Math 428
#30 April 2015

x1=1; #y(0)
x2=0; #y'(0)
x= [x1,x2];
step=0.5; #Step size
boundary = 1;
n=boundary/step;

for i=1:n
  k1=[x2,-sin(x1)];
  k2=[x2+step*k1(1)/2,-sin(x1+step*k1(2)/2)];
  k3=[x2+step*k2(1)/2,-sin(x1+step*k2(2)/2)];
  k4=[x2+step*k3(1),-sin(x1+step*k3(2))];
  x=x+(step/6)*(k1+2*k2+2*k3+k4);
end
x

which gives me x= 1.000000000000000 -0.719893726895089. This feels wrong to me, as x(1)=y(t). Any idea what' s up?

A couple of things. First notice how I edited your post. I didnt change the content, just the appearance of your code. It's now marked as code. You can click on the "+" symbol in the edit menubar, or just do it manually by using code tags.

That's not what I wanted to write about. You created a bug in your code. A number of things contributed to this:

• You used very similar names for things: x1, x2, and x. That's what got you in trouble.
• It looks like you used cut-and-paste and then modified each line.
• You didn't use Matlab to anything close to its full power.
• You didn't use the debugger.

It's been a long time since I've used Matlab, and I don't have it now. So this is untested code.

# First, let's get rid of the potential for the x1 versus x(1) bug.
y0=1; #y(0)
ydot0=0; #y'(0)
x = [y0, ydot0];

step=0.5; #Step size
boundary = 1;
n=boundary/step;

# Define the derivative function.
# Someone wants you to integrate a different function?
# Just change the derivative function.
f = @(x) [x(2), -sin(x(1))];

for i=1:n
  k1=f(x);
  k2=f(x + (step/2)*k1);
  k3=f(x + (step/2)*k2);
  k4=f(x + step*k2);
  x=x+(step/6)*(k1+2*k2+2*k3+k4);
end
x

You could get even more elegant and use the Butcher tableau instead those hard-coded statements inside the loop. Now you can swap out the integration technique for a different Runge-Kutta technique as well as swapping out the function to be integrated.

 

[Haydo]

Yeah, I copy and pasted some old code I wrote for a non-vector-valued Runge-Kutta. Thanks for the tips. What exactly is the debugger? I'm running Octave on Linux; is that still a thing?

 

[D H]

https://www.gnu.org/software/octave/doc/interpreter/Debugging.html