Help me please Plotting values using Secant method

Question

tom balestra 2020-8-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/573769-help-me-please-plotting-values-using-secant-method

评论： Alan Stevens 2020-11-13

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