Help me please Plotting values using Secant method

Hello Everyone,

i wrote a code that describes velocity and position of ball movement on x,y axis and i also have the angular velocity of the ball.

in that code i had to use Runnge Kutta 4-5 order and for the plot, i have to find the time tf with minimum tolerance of 1e-3 while y=0, using the Secant Method.

but i dont know how to start writing the secant method for this one and i hope someone could help me.

clc;
close all;
clear;

g=9.81;
Cd=0.1;
BtI=1;
S=0.0025;
A=0.0415;
p=1.2;
m=0.422;
t0=0;
tf=3;
v=20;
w0=20*pi;
Phi=[10,20,30,45,60];
Tol=1e-3; 

Xdot=@(X2,Y2,T2) (1/m)*(-(S)*T2*Y2-0.5*(Cd)*p*A*X2*(X2^2+Y2^2)^0.5);
Ydot=@(X2,Y2,T2) (1/m)*(-m*g+S*T2*X2-0.5*Cd*p*A*Y2*(X2^2+Y2^2)^0.5);
Tdot=@(X2,Y2,T2) -BtI*T2;
n=1000;
T=linspace(t0,tf,n+1);
h=(tf-t0)/n;

for k=1:length(Phi)
X=zeros(n+1,3);
XD=zeros(n+1,3);
X(1,1:3)=[0,0,0];
XD(1,1:3)=[v*cosd(Phi(k)),v*sind(Phi(k)),w0];


%RK-4
for i=1:n                       
    X2=XD(i,1);
    Y2=XD(i,2);
    T2=XD(i,3);
    
f1=Xdot(X2,Y2,T2);F1=X2;
k1=Ydot(X2,Y2,T2);K1=Y2;
g1=Tdot(X2,Y2,T2);G1=T2;

    X2=XD(i,1)+(h/2)*f1;
    Y2=XD(i,2)+(h/2)*k1;
    T2=XD(i,3)+(h/2)*g1;
    
f2=Xdot(X2,Y2,T2);F2=X2;
k2=Ydot(X2,Y2,T2);K2=Y2;
g2=Tdot(X2,Y2,T2);G2=T2;

    X2=XD(i,1)+(h/2)*f2;
    Y2=XD(i,2)+(h/2)*k2;
    T2=XD(i,3)+(h/2)*g2;
    
f3=Xdot(X2,Y2,T2);F3=X2;
k3=Ydot(X2,Y2,T2);K3=Y2;
g3=Tdot(X2,Y2,T2);G3=T2;

    X2=XD(i,1)+h*f3;
    Y2=XD(i,2)+h*k3;
    T2=XD(i,3)+h*g3;
    
f4=Xdot(X2,Y2,T2);F4=X2;
k4=Ydot(X2,Y2,T2);K4=Y2;
g4=Tdot(X2,Y2,T2);G4=T2;

X(i+1,1)=X(i,1)+(h/6)*(F1+2*F2+2*F3+F4);
X(i+1,2)=X(i,2)+(h/6)*(K1+2*K2+2*K3+K4);
X(i+1,3)=X(i,3)+(h/6)*(G1+2*G2+2*G3+G4);
XD(i+1,1)=XD(i,1)+(h/6)*(f1+2*f2+2*f3+f4);
XD(i+1,2)=XD(i,2)+(h/6)*(k1+2*k2+2*k3+k4);
XD(i+1,3)=XD(i,3)+(h/6)*(g1+2*g2+2*g3+g4);


% y=100; string///////////////////////////////////
% while abs(y) > Tol
%    fb=
%    fa=
%    c=
%    a=b;
%    b=c;
%    abs=a-b;
% end




figure(1)
hold on


end
plot(X(:,1),X(:,2));
end



title('Vo=20','FontSize',18)
ylabel('y','FontSize',18)
xlabel('x','FontSize',18)
ylim([0 30])
grid on

Its supposed to looked like this

but acctualy its looks like this

Please help me

Answer 1

[Alan Stevens]

Increase tf from 3 to 4 to get all the lines returning to zero.

[tom balestra]

yes you are right but i need to do this approximation using Secant method. i just need some help for writing it.

thanks Alan.

[Alan Stevens]

The secant method is for finding the roots of algebraic equations. Do you mean that you have to find the value of tf for each curve that resuts in the curve just returning to zero?

[tom balestra]

exactly, i have to find the time tf with minimum tolerance of 1e-3.

do you have a clue how can i find it?

realy appreciate your efforts

[Alan Stevens]

Well, the secant method here uses two values, say tf1 and tf2 (guessed initially), then updates the next guess using

tf(i) = tf(i-1) - f(tf(i-1))*(tf(i-1) - tf(i-2))/(f(tf(i-1)) - f(tf(i-2)))

where f is the function value: in this case f is the value of the curve at the end time tf.

This value is then used as a new tf2 and the old tf2 becomes the new tf1.

This is repeated until you have convergence.

I suggest you try this for a single angle initially, before extending it to several angles.

Answer 2

[konda sricharan]

matlab code for f(x)=cos(x)-x*exp(x)

[Alan Stevens]

This can be written as a function in the following way

f = @(x) cos(x) - x.*exp(x);

Answer 3

[konda sricharan]

i dont know