1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

RK4 for a system of 2 second order DEs

  1. Apr 1, 2013 #1
    1. For this problem, we are asked to design a MatLab/Octave code that will calculate and plot the orbit of a sattelite



    2. We are given the following 2nd order equations
    [itex] x''= \frac{-x}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]
    [itex]y''= \frac{-y}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]

    We can convert this into the following system of 4 1st order DEs by substitiuting z=x' and w=y'
    [itex]x'=z[/itex]
    [itex]y'=w[/itex]
    [itex]z'= \frac{-x}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]
    [itex]w'= \frac{-y}{( \sqrt{x^{2}+y^{2}})^{3}}[/itex]

    We use the following initial vallues
    [x0 y0 z0 w0] = [4 0 0 0.5]


    3. The attempt at a solution
    I managed to build a program that works well enough that, after adding annotations, I could hand it in. However, as you can see below, the lines within the for loop could use a little help as far as elegance goes (I've attached the .m file if my copy/paste is too difficult to read). Ideally, I'd like to be able to make those lines a little more compact. Any ideas?

    X0=[4 0 0 0.5];
    t0=0;
    tmax=50;
    h=0.25;
    SOL=[t0 X0;];
    f1=inline('z','z');
    f2=inline('w','w');
    f3=inline('-x/((sqrt(x^2+y^2))^3)','x','y');
    f4=inline('-y/((sqrt(x^2+y^2))^3)','x','y');
    for i = 1:tmax/h
    g1=[f1(X0(3)) f2(X0(4)) f3(X0(1),X0(2)) f4(X0(1),X0(2))];
    g2=[f1(X0(3)+0.5*h*g1(3)) f2(X0(4)+0.5*h*g1(4)) f3(X0(1)+0.5*h*g1(1),X0(2)+0.5*h*g1(2)) f4(X0(1)+0.5*h*g1(1),X0(2)+0.5*h*g1(2))];
    g3=[f1(X0(3)+0.5*h*g2(3)) f2(X0(4)+0.5*h*g2(4)) f3(X0(1)+0.5*h*g2(1),X0(2)+0.5*h*g2(2)) f4(X0(1)+0.5*h*g2(1),X0(2)+0.5*h*g2(2))];
    g4=[f1(X0(3)+h*g3(3)) f2(X0(4)+h*g3(4)) f3(X0(1)+h*g3(1),X0(2)+h*g3(2)) f4(X0(1)+h*g3(1),X0(2)+h*g3(2))];
    X1=X0+(h/6)*(g1+2*g2+2*g3+g4);
    t1=t0+h;
    SOL=[SOL;t1 X1];
    X0=X1;
    t0=t1;
    end
    plot(SOL(:,2),SOL(:,3))
    axis('equal')
     

    Attached Files:

    Last edited: Apr 1, 2013
  2. jcsd
  3. Apr 2, 2013 #2
    In case anyone is interested, I did figure something out. It added a few lines to the code, but it did make the loop a bit prettier.

    %sets initial position and velocity values for x and y
    %(z and w are the velecities of x and y, respectively)
    X0=[4 0 0 0.5];
    %initializes the time variable (x, y, z, and w are functions of time)
    %and sets a maximum t for which to approximate
    t0=0;
    tmax=50;
    %stepsize
    h=0.25;
    %initializes array in which to store approximated valies.
    SOL=[t0 X0];
    %these inline functions are our system of 4 1st order DEs, x', y', z', and w'
    f1=inline('z','z');
    f2=inline('w','w');
    f3=inline('-x/((x^2+y^2)^1.5)','x','y');
    f4=inline('-y/((x^2+y^2)^1.5)','x','y');
    %approximates x, y, z, and w in increments of 0.0025 from t=0 to t=50 using the
    %classical 4th order Runge-Kutta method
    for i = 1:tmax/h
    g1=[f1(X0(3)) f2(X0(4)) f3(X0(1),X0(2)) f4(X0(1),X0(2))];
    k1=X0+0.5*h*g1;
    g2=[f1(k1(3)) f2(k1(4)) f3(k1(1),k1(2)) f4(k1(1),k1(2))];
    k2=X0+0.5*h*g2;
    g3=[f1(k2(3)) f2(k2(4)) f3(k2(1),k2(2)) f4(k2(1),k2(2))];
    k3=X0+h*g3;
    g4=[f1(k3(3)) f2(k3(4)) f3(k3(1),k3(2)) f4(k3(1),k3(2))];
    X1=X0+(h/6)*(g1+2*g2+2*g3+g4); %RK4
    t1=t0+h;
    SOL=[SOL;t1 X1];
    X0=X1;
    t0=t1;​
    end
    %plots our approximation to the position functions, x and y and sets the axes
    %to use equal scales
    plot(SOL(:,2),SOL(:,3))
    axis('equal')
     
    Last edited: Apr 2, 2013
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted