How to use anonymous functions and fsolve to find initial conditions given final conditions (projectile motion)

What I want to do here is rework my initial code to get the initial velocity and angle given a target distance (x,y).

this is based off of a system that uses 4 springs to launch a ball. therefore inital velocity is mainly dependent on the pull length L.

code is as follows:

 %PULL LENGTH
L=0.2;
%LAUNCH ANGLE
theta(1)=45;
q(2)=theta(1);
%INITIAL POSITION OF LAUNCH
x(1)=0;
y(1)=0;
 %INITIAL VELOCITY
g=9.81;
%MASS OF BALL
mb=0.0116;
%MASS TOTAL
mt=0.08193;
%SPRING CONSTANT
k=42; %N*m IN EACH OF 4 SPRINGS
%PRE-LOAD
p=0.0125;
%INITIAL VELOCITY AND COMPONENTS
v(1)=sqrt((4*k*(L^2-p^2)-2*mt*g*L*sin(theta(1)))/mt);
q(1)=v(1);
vx(1)=v(1)*cos(theta(1));
vy(1)=v(1)*sin(theta(1));
%DRAG
d=(0.01588*2);
A=pi()*d^2/4;
rho=1.2;
cd=0.5;
D=1/2*cd*rho*A;
%TRAJECTORY LOOP
del_t=0.002;
for n=1:2000;
t(n)=n*del_t;
v(n)=sqrt(vx(n)^2+vy(n)^2);
theta(n)=atan(vy(n)/(vx(n)));
Fd(n)=D*v(n)^2;
ax(n)=-Fd(n)*cos(theta(n))/mb;
ay(n)=-Fd(n)*sin(theta(n))/mb-g;
x(n+1)=x(n)+vx(n)*del_t;
y(n+1)=y(n)+vy(n)*del_t;
vx(n+1)=vx(n)+ax(n)*del_t;
vy(n+1)=vy(n)+ay(n)*del_t;
F = @LAUNCH;
q = fsolve(F,q(1));
end
plot(x,y);

with this function in a separate file

 function [F] = LAUNCH(q);
    F(1)= q(1)-sqrt((4*k*(L^2-p^2)+2*mt*g*L*sin(theta(1)))/mt);
    F(2)= mt*g*L*sin(q(2));
end

Im assuming here that I cant use x(1) and x(2) in the function because they are already variables in the original code.

What is have now is not functioning properly and q returns v1 correctly and also the angle I input. How do I make this loop so that it gives me v1 and theta such that it will hit a target?

[Tamir Suliman]

you could provide more explanantion on what equation you are confused

i think you could use a function of x,y and inline to sub the values of x ,y

example

 F =@(x) (-2*x.^6)-(1.6*x.^4)+(12*x)+1;

[Nathan Stone]

What would this look like in practice? Im trying to figure out how anonymous functions and fsolve work but Im having difficulty.

so far I have

 Function F = LAUNCH(x)
    F(1)=v(1)-sqrt((4*k*(L^2-p^2)+2*mt*g*L*sin(theta(1)))/mt)
    F(2)=

now I need another function relating to theta correct? Or am I going about this the wrong way?

[Daniel kiracofe]

What I was suggesting was to take more or less your entire code, wrap that into a function, and then solve it. i.e.

function xf = launch( L, theta )
 %PULL LENGTH
%  L=0.2; %now an input from the function
%LAUNCH ANGLE
%  theta(1)=45;  %now an input from the function
q(2)=theta(1);
%INITIAL POSITION OF LAUNCH
x(1)=0;
y(1)=0;
 %INITIAL VELOCITY
g=9.81;
%MASS OF BALL
mb=0.0116;
%MASS TOTAL
mt=0.08193;
%SPRING CONSTANT
k=42; %N*m IN EACH OF 4 SPRINGS
%PRE-LOAD
p=0.0125;
%INITIAL VELOCITY AND COMPONENTS
v(1)=sqrt((4*k*(L^2-p^2)-2*mt*g*L*sin(theta(1)))/mt);
q(1)=v(1);
vx(1)=v(1)*cos(theta(1));
vy(1)=v(1)*sin(theta(1));
%DRAG
d=(0.01588*2);
A=pi()*d^2/4;
rho=1.2;
cd=0.5;
D=1/2*cd*rho*A;
%TRAJECTORY LOOP
del_t=0.002;
for n=1:2000;
t(n)=n*del_t;
v(n)=sqrt(vx(n)^2+vy(n)^2);
theta(n)=atan(vy(n)/(vx(n)));
Fd(n)=D*v(n)^2;
ax(n)=-Fd(n)*cos(theta(n))/mb;
ay(n)=-Fd(n)*sin(theta(n))/mb-g;
x(n+1)=x(n)+vx(n)*del_t;
y(n+1)=y(n)+vy(n)*del_t;
vx(n+1)=vx(n)+ax(n)*del_t;
vy(n+1)=vy(n)+ay(n)*del_t;
end
vf = vy(end)

Then for a given pull length and angle, you have a function that gives the final position. You then use fsolve to solve for pull length and angle the correspond to the desired final position. Computationally might not be the quickest way to do it, but it's straightfoward (once you know how fsolve works). Walter's suggestion to integrate the ODE backwards in time is almost certainly computationally quicker. That one requires more math thinking and less programming thinking. My idea requires less math but more programming thinking.

Answer 1

[Daniel kiracofe]

So it looks like you are integrating an ODE numerically. Starting from a set initial conditions, you are doing a simple time integration. If I understand the problem correctly, what you want to is to solve the inverse problem: give a set of final conditions, find a set of initial conditions that get there.

Few ways that you could do it.

1) Instead of integrating numerically, solve the ODE analytically. This function will explicitly contain the initial conditions, and so can be easily inverted. If there was no drag, this would be pretty easy: x(t) = x(0) + vx(0) * t, y(t) = y(0) + vy(0) * t + (1/2) g * t^2 (or something reasonably close to that). Then just substitute in the final conditions that you want, and solve for vx(0) and vy(0). The presence of drag makes the function more complicated, but if you are clever you can figure it out. Symbolic math toolbox, or perhaps a program like Maple or Mathematica might be able to help.

2) If you want to do it numerically, one option is to basically be doing is a guess and check strategy. i.e. guess an initial condition, run through your code, see how close it is, try another guess, repeat. Matlab can automate this for you with function like lsqnonlin() or fsolve(). As Tamir suggested, create a function from your code. i.e. the function would take as an input the initial conditions and give as an output the final conditions. then pass that entire function to fsolve and it will solve it for you.

[Walter Roberson]

To determine the set of initial conditions that got you to a final situation, you can run an ode backwards in time.

[Nathan Stone]

Youre right, Im looking to take final values of distance (x and y) and solve for a set of initial conditions, in this case launch angle and velocity.

I believe that in terms of simplicity your second option is probably the better one. I knew there was a way to do a "guess and check" but Im unsure how to go about this exactly, im pretty new to matlab.

So for the function I would input something like

 [Final_x_y] = function(theta,v(n))

then use fsolve to find theta and v(n)?