BVP_projectile_with_drag_sol_k

% 抛体运动——用bvp求解阻力系数
g=10;k=1;theta=pi/4;n=2;v=10;
solinit = bvpinit(linspace(0,0.734846711099164,1e4), ...
                  @(t)mat4init(t,g,v,theta),k);
opts_jac = bvpset('FJacobian',@(t,u,k)jac(t,u,k,n), ...
             'BCJacobian',@(u0,ut,k)bc_ab(u0,ut,k), ...
             'RelTol',0.1,'AbsTol',0.1,'Stats','on');
opts = bvpset('RelTol',0.1,'AbsTol',0.1,'Stats','on');
tic
sol = bvp4c(  @(t,u,k)mat4ode(t,u,k,n,g), ...
            @(u0,ut,k)mat4bc(u0,ut,k,v,theta), solinit,opts_jac);
toc
% 输出验证
fprintf('\n阻力系数k=：%7.3f.\n',sol.parameters)
t=linspace(0,0.7348,1e3)';
ut = deval(sol,t)';
x=ut(:,1);y=ut(:,2);
plot(x,y)
% 微分方程
function du = mat4ode(~,u,k,n,g)
vx = u(3); vy = u(4);
v  = sqrt(vx^2+vy^2);
dx  = vx;
dy  = vy;
dvx =   -k*v^(n-1)*vx;
dvy = -g-k*v^(n-1)*vy;
du  = [dx dy dvx dvy]';
end
% FJacobian
function [J,du_k]= jac(~,u,k,n)
vx = u(3); vy = u(4);
v  = sqrt(vx^2+vy^2);
dvx_vx = -k*v^(n-1)-k*(n-1)*v^(n-3)*vx*vx;
dvx_vy = -k*(n-1)*v^(n-3)*vx*vy;
dvy_vx = -k*(n-1)*v^(n-3)*vx*vy;
dvy_vy = -k*v^(n-1)-k*(n-1)*v^(n-3)*vy*vy;
J = [
    0 0 1 0;
    0 0 0 1;
    0 0 dvx_vx dvx_vy;
    0 0 dvy_vx dvy_vy];
du_k=[0 0 -v^(n-1)*vx -v^(n-1)*vy]';
end
% 边界条件
function bc = mat4bc(u0,ut,~,v,theta)
x0=u0(1);y0=u0(2);vx0=u0(3);vy0=u0(4);
xt=ut(1);yt=ut(2);vxt=ut(3);vyt=ut(4);
ua=[x0 y0 vx0-v*cos(theta) vy0-v*sin(theta)]';
bc = [ua(1:4);yt];
end
% BCJacobian
function [a,b,k]=bc_ab(~,~,~)
a=[ 1 0 0 0 0
    0 1 0 0 0
    0 0 1 0 0
    0 0 0 1 0
    0 0 0 0 0 ];
b=[ 0 0 0 0 0
    0 0 0 0 0
    0 0 0 0 0
    0 0 0 0 0
    0 1 0 0 0 ];
k=[ 0 0 0 0 0 ]';
end
% 初值估计
function u_t = mat4init(t,g,v,theta)
vxt=v*cos(theta);vyt=v*sin(theta);
u_t=[vxt*t vyt*t-0.5*g*t^2 vxt vyt]';
end