Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab programming using shooting method, Euler and Runge Ketta methods

  1. Nov 29, 2011 #1
    I need help to solve this coursework:

    MATLAB PROGRAMMING COURSEWORK
    OBJECTIVES:
     Learn to solve engineering problems using MATLAB
     Write Euler and Runge-Kutta initial-value ODE solvers
     Write a Shooting Method boundary-value ODE solver
     Investigate the properties of the solvers
     Summarise your work in a short report
    2. Assignment Part 2 – Writing a shooting method BVP ODE solver
    A test of the Beagle 3 Mars rover is being conducted (on Earth). The rover is
    designed to enter the Mars atmosphere, deploy two parachutes, inflate
    airbags and bounce on impact before coming to rest. For this test, the rover
    will be dropped from beneath an aircraft and its descent will be filmed, so the
    release and landing points must be precise. The rover will be released 10km
    above the surface, and the target forward distance from the release point is
    1km. Your job is to calculate the required velocity of the aircraft at the point of
    release, and you will achieve this using the Shooting Method. Further
    information is provided below:
    Rover
    Mass [kg] 50
    Frontal area [m2] 0.5
    Drag coefficient 0.7
    Pilot parachute
    Release time 60
    Area [m2] 1.0
    Drag coefficient 1.0
    Main parachute
    Release time 120
    Area [m2] 5
    Drag coefficient 1.0
    Airbags
    Inflation height [m] 200
    Frontal area [m2] 3.5
    Drag coefficient 0.8
    Stiffness [N/m] 10000
    Damping [Ns/m] 300
    Other Air density at sea level [kg/m3] 1.207
    Suggested plan of attack:
    1. Work out the ODE for this system and convert it to state-space form.
    2. Modify your Dz.m function from part 1 to calculate the state derivatives.
    3. Modify your odeSolver code to compute and plot the trajectory of the
    rover, given initial values for the states. You should mark the points of
    parachute release. Run this to make sure it is working before moving
    on to the shooting method.
    4. Create a new function to use the shooting method to calculate the required
    aircraft velocity to hit the target distance. This function should repeatedly
    call all of the modified functions you have just created. Check the notes
    to see how the shooting method works. Think carefully about how
    your new function will work before starting – sketch it out.
    Just like any real Engineering problem, you will need to make (and justify)
    simplifying assumptions to solve the problem.
    Once you have solved the problem – well done! However, this is the
    minimum requirement to pass the coursework. Successful completion of
    this work, combined with neat code and a good report could get you into 2:1
    territory. Access to higher marks is gained by taking the problem further. This
    exercise is open-ended – what you choose to explore is entirely up to you.
    Here are some suggestions to get you started:
     Think about the simplifying assumptions you have made in arriving at
    your ODE. How can you make the solution more realistic? For
    example, the air density changes with height above the ground. How
    does this affect the problem?
     What about making the target distance a user input?
     You may have discovered that it is easy to make your code fail to reach
    a solution. Investigate what causes it to fail and find ways to make it
    more robust.
     How could you make your code more accurate or efficient?
    Write a report of maximum length 4 pages containing:
     The derivation of your ODE
     A graph of the trajectory, showing the point of parachute deployment
     The impact speed
     Discussion of simplifying assumptions
     Discussion of the further work you have completed
     Your commented code included in an appendix (not included in the 4 pages limit)
     
  2. jcsd
  3. Nov 29, 2011 #2

    hunt_mat

    User Avatar
    Homework Helper

    We won't do your homework for you but I will suggest a few good places to start.

    Write down the differential equation that you need to solve is a good place to start. What is your differential equation.
     
  4. Nov 30, 2011 #3
    Could you please check my last ODE hereafter?

    function [dz] = Dz2(t, z)
    %initial conditions
    persistent hasbounced;
    g=9.81; %gravity constant
    M=50; %Mass of rover
    R=300; %Damper coefficient
    k=10000; %spring stiffness
    V=((z(2))^2)+((z(4))^2); %Velocity
    theta=atan2((z(4)),(z(2)));
    %density
    d=Density(z(3));


    if t<60
    hasbounced=0;
    a=0.5; %a=area
    c=0.7; %c=drag coefficient
    elseif t<120
    a=1.0;
    c=1.0;
    else
    a=5.0;
    c=1.0;
    end
    %conditions for if hasbounded is true
    if hasbounced
    a=3.5;
    c=0.8;
    end
    Fd=(1/2)*d*V*c*a;

    dz(1) = z(2);
    dz(2) = -(Fd*(cos(theta)))/M;
    dz(3) = z(4);

    if z(3)>1.06 %the radius of the airbag
    dz(4) = -g-(Fd*(sin(theta)))/M;
    else
    dz(4) = -g+(((k*(0.66-(z(3)-0.4)))-(R*z(4))))/M; %0.4 is the radius of the rover
    hasbounced=1;
    %hasbounced is true only after programme has run new dz(4) process for first times
    end

    dz=dz';

    I don't know how to add it the friction when the Rover hits the ground.

    Can you help me as well on this topic?
     
  5. Dec 6, 2011 #4
    Could you please check my last ODE hereafter?

    function [dz] = Dz2(t, z)
    %initial conditions
    persistent hasbounced;
    g=9.81; %gravity constant
    M=50; %Mass of rover
    R=300; %Damper coefficient
    k=10000; %spring stiffness
    V=((z(2))^2)+((z(4))^2); %Velocity
    theta=atan2((z(4)),(z(2)));
    %density
    d=Density(z(3));


    if t<60
    hasbounced=0;
    a=0.5; %a=area
    c=0.7; %c=drag coefficient
    elseif t<120
    a=1.0;
    c=1.0;
    else
    a=5.0;
    c=1.0;
    end
    %conditions for if hasbounded is true
    if hasbounced
    a=3.5;
    c=0.8;
    end
    Fd=(1/2)*d*V*c*a;

    dz(1) = z(2);
    dz(2) = -(Fd*(cos(theta)))/M;
    dz(3) = z(4);

    if z(3)>1.06 %the radius of the airbag
    dz(4) = -g-(Fd*(sin(theta)))/M;
    else
    dz(4) = -g+(((k*(0.66-(z(3)-0.4)))-(R*z(4))))/M; %0.4 is the radius of the rover
    hasbounced=1;
    %hasbounced is true only after programme has run new dz(4) process for first times
    end

    dz=dz';

    I don't know how to add it the friction when the Rover hits the ground.

    Can you help me as well on this topic?
     
  6. Dec 6, 2011 #5
    http://www.bath.ac.uk/mech-eng/people/profiles/hillis/hillis.jpg [Broken]

    This guy could help.
     
    Last edited by a moderator: May 5, 2017
  7. Dec 9, 2011 #6
    hello.... i need help :P
    this is the "code"

    function [x,w1,w2]=disnolin(f,a,b,d,c,n,tol,m)
    %paso1
    w1=zeros(1,n+1);
    w2=zeros(1,n+1);
    x=zeros(1,n+1);
    f1=inline(diff(f,'w'),x,w,z);
    f2=inline(diff(f,'z'),x,w,z);
    h=(b-a)/n;
    k=1;
    TK=(c-d)/(b-a);
    %paso2
    while k<=m
    %paso3
    w1(1)=d;
    w2(1)=TK;
    u(1)=0;
    u(2)=1;
    %paso4
    for i=2:n
    %paso5
    x=a+(i-2)/h;
    %paso6
    k11=h*w2(i-1);
    k12=h*f(x,w1(i-1),w2(i-1));
    k21=h*(w2(i-1)+0.5*k12);
    k22=h*f(x+0.5*h,w1(i-1)+0.5*k11,w2(i-1)+0.5*k12);
    k31=h*(w2(i-1)+0.5*k22);
    k32=h*f(x+0.5*h,w1(i-1)+0.5*k21,w2(i-1)+0.5*k22);
    k41=h*(w2(i-1)+k32);
    k42=h*f(x+h,w1(i-1)+k31,w2(i-1)+k32);
    w1(i)=w1(i-1)+(k11+2*k21+2*k31+k41)/6;
    w2(i)=w2(i-1)+(k12+2*k22+2*k32+k42)/6;
    l11=h*u(2);
    l12=h*(f1(x,w1(i-1),w2(i-1))*u(1)+f2(x,w1(i-1),w2(i-1))*u(2));
    l21=h*(u(2)+0.5*l12);
    l22=h*(f1(x+0.5*h,w1(i-1),w2(i-1))*(u(1)+0.5*l11)+f2(x+0.5*h,w1(i-1),w2(i-1))*(u(2)+0.5*l12));
    l31=h*(u(2)+0.5*l22);
    l32=h*(f1(x+0.5*h,w1(i-1),w2(i-1))*(u(1)+0.5*l21)+f2(x+0.5*h,w1(i-1),w2(i-1))*(u(2)+0.5*l22));
    l41=h*(u(2)+l32);
    l42=h*(f1(x+h,w1(i-1),w2(i-1))*(u(1)+l31)+f2(x+h,w1(i-1),w2(i-1))*(u(2)+l32));
    u(1)=u(1)+(l11+2*l21+2*l31+l41)/6;
    u(2)=u(2)+(l12+2*l22+2*l32+l42)/6;
    %paso7
    if abs(w1(n)-c)<=tol
    %paso8
    for j=0:n
    x=a+j*h;
    %[x' w1(i)' w2(i)'];
    %paso9
    end
    end
    end
    %paso10
    TK=TK-(w1(i)-c)/u(1);
    k=k+1;
    end
    %paso11
    fprint('numero maximo de ieraciones excedido');




    function y=f(x,w,z)
    y=(32+2*x.^3-w.*z)./8;


    thx!!!
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Matlab programming using shooting method, Euler and Runge Ketta methods
  1. Euler method (Replies: 12)

  2. Euler's method (Replies: 5)

  3. Euler's Method (Replies: 3)

Loading...