Heun's method program code

Hello,

I am trying to program a script to solve a second order ODE using the Heun's method as required for a project of mine. I cannot use ODE45 or any of the like. Here is the code I have written so far:

 function project2()

 h = 1/32; %STEP SIZE
z0 = [0, 2.4]; %INITIAL Y VALUE

 t=0:h:20;
 t(1)=0; 
 z(1)=0;
 n=1;

 while t < 6 %TIME INTERVAL
    t(n+1) = t(n) + h;
    k1 = f(t(n), z(n));
    k2 = f(t(n)+h, z(n)+k1*h);
    z(n+1) = z(n) + (h/2)*(k1+k2);
    n = n+1;
end

 plot(t,z,'g-','LineWidth',1)

 function dzdt = f(t,z)
w = 2;
G = 9.81;
z1=z(1); % get z1
z2=z(2); % get z2
dzdt = [z2 ; -G*sin(2*t)+w^2*z1;];

I run this and get no errors, but the lot comes up completely empty. Please help!

Answer 1

[Geoff Hayes]

Sean - the problem is with your while loop condition

 while t < 6

At this point, t is an array of 641 elements because of the line

 t=0:h:20;

It is unclear why you populate this array as such (with the step size of h) and then reset the first element to 0

 t(1)=0;

and then reset every element of t on the first line of the while loop body. You don't need to do this if t has already been initialized correctly.

Instead of the above code, try using a for loop with n incrementing on each iteration of the loop. Perhaps

   function project2()
   h = 1/32; %STEP SIZE
   z0 = [0, 2.4]; %INITIAL Y VALUE
   t=0:h:20;
   z(1)=0;
   for n=1:length(t)-1
      k1 = f(t(n), z(n));
      k2 = f(t(n)+h, z(n)+k1*h);
      z(n+1) = z(n) + (h/2)*(k1+k2);
  end

The above code will now call your f function, but there is a bug with that too

function dzdt = f(t,z)
w = 2;
G = 9.81;
z1=z(1); % get z1
z2=z(2); % get z2
dzdt = [z2 ; -G*sin(2*t)+w^2*z1;];

The code for this function assumes that z is two dimensional...but you only supply a scalar when calling it from within your loop

k1 = f(t(n), z(n));
k2 = f(t(n)+h, z(n)+k1*h);

What should you be supplying to this function f? z(n-1:n) to get the two elements? Please clarify.

[Sean Malinowski]

Geoff, thank you for the help so far! I figured that my function was not correct. Honestly, I do not know what it should look like. I have practically no MATLAB experience or training and was all of a sudden assigned to create this program. What I have so far is from hours upon hours of scrounging the internet looking for answers. I have read that the function must be put into a column vector and I know that the notation for z(n) should be different, but I do not know what it should be.

[Geoff Hayes]

Sean - can you provide some details on your f function? What are you basing it on?

[Sean Malinowski]

Geoff,

The function is supposed to represent the second order differential equation:

 r''=-Gsin(wt)+w^2r

"r" can be replaced with y respectively, this is just the dependent variable that was noted in my project. Initial conditions are:

 r(0)=0
 r'(0)=v0

v0 is the initial velocity and will be changing from 2.4, 2.45 and 2.5. I attempted to convert this into a system of first order differential equations, hence the code:

 z1=z(1); % get z1
 z2=z(2); % get z2
 dzdt = [z2 ; -G*sin(2*t)+w^2*z1;];

But, I do not know if I did this correctly as I have only written code for single first order ODEs.

Answer 2

[ayman alarousi]

i want mathlab codes for second order ODE

midpoint, Runge-Kutta method

[Torsten]

https://de.mathworks.com/matlabcentral/answers/159294-midpoint-numerical-integration-without-a-built-in-function