2nd Order Adams-Bashforth/Runge Kutta

[CMJ96]

Homework Statement

s=0.140406704
2. Relevant equation

The Attempt at a Solution

So I converted the ODE into the following two equations
$$\frac{dx_1}{dt}=x_2$$
$$\frac{dx_2}{dt}=x_1- \alpha sin(x_2)$$

I have done the following with the program so far, I came to a halt because I am not sure how to define k1 and k2 with functions of 3 variables, also I find it a little confusing that there is a time-step despite there being no t's in the functions?
program ab2
implicit none
integer :: i,n
real :: t,x1,x2,func1,func2,h, t0, x10,x20,s,k11,k12,k22,k21
a=0.5
h=0.01
t0=0
t=t0
s=0.140406704
x10=s
x20=-s
x1=x10
x2=x20

end program

function func1(x1,x2,t)
real :: func1,x1,t,x2
func1= x2
end function func1

function func2(x1,x2,t)
real :: func2,x1,x2,t,a
func2= x1-a*sin(x2)
end function func2

 

[BvU]

Please enclose source code in [code=fortran]...[/code] tags:

program ab2
implicit none
integer :: i,n
real :: t,x1,x2,func1,func2,h, t0, x10,x20,s,k11,k12,k22,k21
a=0.5
h=0.01
t0=0
t=t0
s=0.140406704
x10=s
x20=-s
x1=x10
x2=x20

end program

function func1(x1,x2,t)
real :: func1,x1,t,x2
func1= x2
end function func1

function func2(x1,x2,t)
real :: func2,x1,x2,t,a
func2= x1-a*sin(x2)
end function func2

there is a time-step despite there being no t's in the functions

so $f_1 = f_1(x_1, x_2) = f_1(x_2)$ and $f_2=f_2(x_1,x_2)$ and the time dependence comes in when you take the step. What I miss in your account is the realization that there is one $f$ with 2 components: the first is the derivative of $x_1$ and the second is the derivative of $x_2$. It's ok to call the first $f_2$ and the second $f_2$.

Perhaps you can first work this out with a first order integrator (Euler = first order Runge Kutta) and then expand to second order ?[edit] oh, and: how wdoes func2 know about a ?

 

[CMJ96]

So would it be wise to remove the t variables from the functions, and the x1 variable from function 1?
Also, I've edited the code to define a=0.5 inside the function.

It's ok to call the first f2f2f_2 and the second f2f2f_2.

I'm a bit confused by this, why is it okay to do that?

 

[BvU]

You have a system of two first order differential equations (one for $u = x_1$ and one for $v = x_2$ *) with two initial conditions.
For the step you apply $f_1$ to $x_1$ and $f_2$ to $x_2$

*) or you can consider $u$ as a vector with two components, in which case $k_1$ and $k_2$ are also vectors with two components.

 

[CMJ96]

You have a system of two first order differential equations (one for $u = x_1$ and one for $v = x_2$ *) with two initial conditions.
For the step you apply $f_1$ to $x_1$ and $f_2$ to $x_2$

Does this mean each function will be $f_1(x_1,x_2)$ and $f_2(x_1,x_2)$ and then can I apply Runge-Kutta to those functions, where $x_2$ is equivalent to $t$ in the equation I gave in my initial post for $f_1$ and vise versa for $f_2$?
And because $f_1(x_2)$ the Runge-Kutta method would simplify down quite a lot?

 

[BvU]

We get confused by the names and the subscripts. $x_2$ is definitely not equivalent to $t$. You had a system of two differential equations $$
\begin{align*}
\frac{dx_1}{dt}&=x_2\\ \\
\frac{dx_2}{dt}&=x_1- \alpha \sin(x_2)
\end{align*}\\ $$ and you are going to integrate both $x_1$ and $x_2$ as functions of $t$. You called the subroutine that calculates $dx_1\over dt$ func1 and that is a function of $x_2$ only. It will help you to calculate $k_1$ for $x_1$ at $t=0$. You can not calculate $k_2$ yet, because you do not have the $x_2$ for $t=h$.

So you first need the subroutine that calculates $dx_2\over dt$, named func2 and that is a function of $x_1$ and $x_2$ only.

Try to do one integration step by hand: write it out for $x_1$ and write it out for $x_2$. Then the pattern emerges, I hope...

 

[CMJ96]

We get confused by the names and the subscripts. $x_2$ is definitely not equivalent to $t$. You had a system of two differential equations $$
\begin{align*}
\frac{dx_1}{dt}&=x_2\\ \\
\frac{dx_2}{dt}&=x_1- \alpha \sin(x_2)
\end{align*}\\ $$ and you are going to integrate both $x_1$ and $x_2$ as functions of $t$. You called the subroutine that calculates $dx_1\over dt$ func1 and that is a function of $x_2$ only. It will help you to calculate $k_1$ for $x_1$ at $t=0$. You can not calculate $k_2$ yet, because you do not have the $x_2$ for $t=h$.

So you first need the subroutine that calculates $dx_2\over dt$, named func2 and that is a function of $x_1$ and $x_2$ only.

Try to do one integration step by hand: write it out for $x_1$ and write it out for $x_2$. Then the pattern emerges, I hope...

Okay so I have just worked it out on paper and I think I understand this more now, we couldn't find $k_2$ for the first ODE, because it depends on $x_2$, so we have to go through the second ODE to find $x_2$ at $t+h$, using Runge-Kutta on the second ODE I found it to be -0.140214204 at $t+h$, then I subbed that into the equation for $k_2$ in the first ODE. This got values of $x_1$ and $x_2$ at $t+h$ of 0.140266395 and -0.140214204. This seems to be a reasonable estimate for those values?

 

[BvU]

I don't have a feeling for values. Much better to write it out using symbols/variable names. For the values, the computer does the work. For us, it's about the algorithm: steps to come to a correct program, working as designed.

worked it out on paper

that way I can't comment

[edit]in the sense: the values are not correct, but I don't know what went wrong.

 

[CMJ96]

Is this more like it? This is to calculate $x_{1,1}$ from $x_{0,1}$ and $x_{1,2}$ from $x_{0,2}$.
$$k_{1,1}=h f_1(x_{0,2})$$
$$k_{1,2}=h f_2(x_{0,1},x_{0,2})$$
$$k_{2,1}=h f_1(x_{0,2} + k_{1,2})$$
$$k_{2,2}=h f_2(x_{0,1}+k_{1,1}, x_{0,2} + k_{1,2} ) $$
$$x_{1,1}=x_{0,1} + \frac{1}{2} \left( k_{1,1} + k_{2,1} \right)$$
$$x_{1,2}=x_{0,2} + \frac{1}{2} \left( k_{1,2} +k_{2,2} \right) $$

 

[BvU]

Looks good. The exercise wants to switch over to Adam-Bash from then on.
Does your fortran knowledge include arrays ?

 

[CMJ96]

Looks good. The exercise wants to switch over to Adam-Bash from then on.
Does your fortran knowledge include arrays ?

I have a rough idea on how to use arrays in fortran, I'm not completely comfortable but I can give it a go , for previous exercises involving Adams-Bashforth I've used holding variables and generated text documents of data using the append feature

 

[CMJ96]

This is my progress so far in implementing this into a program

program coupled
implicit none
integer :: i,n
real :: t,func1,func2,h, t0,a,s,k11,k12,k22,k21
real,dimension (25000):: x1,x2
a=0.5
h=0.001
t0=0
t=t0
s=0.140406704
x1(1)=s
x2(1)=-s

n=25000

k11=h*func1(x2(1))

k12=h*func2(x1(1),x2(1))

k21=h*func1(x2(1)+k12)

k22=h*func2(x1(1)+k11,x2(1)+k12)

x1(2)=x1(1)+(1.0/2.0)*(k11+k21)
x2(2)=x2(1)+(1.0/2.0)*(k12+k22)

do i=2,n
x2(i+1)=x2(i)+ (1.0/2.0)*h*(3.0*func2(x1(i),x2(i))-func2(x1(i-1),x2(i-1)))
x1(i+1)=x1(i)+ (1.0/2.0)*h*(3.0*func1(x2(i)))
end do
print*,x2(15000)
print*,x1(15000)
end program

function func1(x2)
real :: func1,x2
func1= x2
end function func1

function func2(x1,x2)
real :: func2,x1,x2,a
a=0.5
func2= x1-a*sin(x2)
end function func2

does this look along the right lines?

 

[BvU]

Well, does it give a result as expected ? I don't immediately spot an error, but that is no guarantee whatsoever

Now about the programming aspects: I suggested something in the direction of using arrays and you did. Just not in the dimension I intended, but in the time direction. Is there a good reason to keep all 25000 $x_1$ and $x_2$ ? In post #1 it seems printing every 50th result is all that's asked.

The 'direction' I intended is the numbering of the equations: you can set up an integration with only one subroutine for the derviative that returns the derivative vector -- thus getting a program that can easily be modified to deal with a system of any number of equations. And you don't drown in the subscripts: there is one $\vec x$ (or $\vec u$ as it's called in post #1) , one $\vec f$ (that is called twice), one $\vec k_1$ and one $\vec k_2$ for RK2.

For AB $\vec f$ is called only once per step, but you'll need to store the preceding $\vec f$. (So AB -- which I didn't know until now -- is cheaper in computation effort, at the cost of requiring a little more storage).

 

[CMJ96]

every 50th result is all that's asked

Could I implement an if statement like the following?

if ( mod(i,50)==0) then

open(unit=1,file='C:\Users\Callum\Documents\Numerical Methods\coupledx1.txt')
write(1,*)x1
close(1)
open(unit=3,file='C:\Users\Callum\Documents\Numerical Methods\coupledx2.txt')
write(3,*)x2
close(3)
end if

 

[BvU]

1) There is a way to find out ... do it !
2) Opening is usually done only once at the beginning. Idem closing at the end!
3) it would be more sensible to write t, x1 and x2 on one line for every desired t in one file.
4) is t = 0 being written out this way ?