I'm trying to do one of MIT's Computational Methods in Aerospace Engineering projects. I've written the code for Project #1 part 1, and I'm getting a very reasonable solution using the euler, midpoint, heun, and classic RK methods. However, when I run an error analysis on these data sets compared to a really small dt classic RK run, the orders of the errors don't match up. Starting with midpoint, i'm getting a max error of around 0.055 for my euler, midpoint, and huen methods. This is telling me that my midpoint and heun methods are somehow not increasing in order from the euler method, and I can't for the life of me figure out what is going on. Here's the code for my midpoint method:
t_val=zeros(1,n); %initial values
h_val=zeros(1,n);
H_val=zeros(1,n);
A_val=zeros(1,n);
a_val=zeros(1,n);
a_val(1)=a0; %0<a<0.08 radians
for i=2:n
t_val(i)=t_val(i-1)+dt;
%{
ka_2=(dt/2)*A_val(i-1);
kh_2=(dt/2)*H_val(i-1);
kA_1=dt*A1(a_val(i-1),A_val(i-1),h_val(i-1),H_val(i-1));
kH_1=dt*H1(A_val(i-1),a_val(i-1),h_val(i-1),H_val(i-1));
kA_2=dt*A1(a_val(i-1)+ka_2,A_val(i-1)+(kA_1/2),h_val(i-1)+kh_2,H_val(i-1)+(kH_1/2));
kH_2=dt*H1(A_val(i-1)+(kA_1/2),a_val(i-1)+ka_2,h_val(i-1)+kh_2,H_val(i-1)+(kH_1/2));
H_val(i)=H_val(i-1)+kH_2;
A_val(i)=A_val(i-1)+kA_2;
a_val(i)=a_val(i-1)+2*ka_2;
h_val(i)=h_val(i-1)+2*kh_2;
%}
t_val is the time array
h_val is the plunge data point array (trying to solve for h(t))
a_val is the pitch data point array (trying to solve for a(t))
H_val is the first derivative of h(t) data point array
A_val is the first derivative of a(t) data point array
A1(a,A,h,H) and H1(A,a,h,H) are the second derivatives of a(t) and h(t) respectively. They are functions of four variables, and they are found by solving the original system outlined in the project definition for the second derivatives.
The "k" notation is based on an example of the midpoint method I found in a textbook. These "k" variables are just temporary stepping stones, building the real solution at the end.
In one of the solutions for the project, the writer provides a discretization for the midpoint method as follows:
I tried to re-write my code based on this, but it confused me, because it uses three instances of the iteration, n. It uses n-1, n, and n+1. How can I write a code that can use all three of these indexes? What is wrong with my code in the first place? General information always welcome.
Your first method seems to be tailored for a system $\dot x=X$, $\dot X=F(x,X)$. Let's shorten that to $\dot y=f(y)$. The method you apply is the explicit midpoint method \begin{align} k_1&=f(y)\,dt, \\ k_2&=f(y+\tfrac12k_1)\, dt, \\ \hline y_{+1}&=y+k_2 . \end{align} Detailed to the components $x$ and $X$ this gives \begin{align} k_{1x}&=X\,dt,& k_{1X}&=F(x,X)\,dt,\\ k_{2x}&=(X+\tfrac12k_{1X})\,dt,& k_{2X}&=F(x+\tfrac12k_{1x},X+\tfrac12k_{1X})\,dt,\\ \hline x_{+1}&=x+(X+\tfrac12k_{1X})\,dt,&X_{+1}&=X+k_{2X}. \end{align} This reconstruction of the method exhibits that your computation of the final updates of $a$ and $h$ is incomplete, using only an order 1 correct increment which renders the whole method as implemented only order 1, even if the theoretical method is order 2.