Matlab programming using shooting method, Euler and Runge Ketta methods

[haiisrael]

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)

 

[hunt_mat]

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.

 

[haiisrael]

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?

 

[haiisrael]

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?

 

[RungeKutta]

http://www.bath.ac.uk/mech-eng/people/profiles/hillis/hillis.jpg

This guy could help.

 

[losky88]

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!