# BVP_projectile_with_drag_sol_k

Mar 29th, 2022
887
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. % 抛体运动——用bvp求解阻力系数
2. g=10;k=1;theta=pi/4;n=2;v=10;
3. solinit = bvpinit(linspace(0,0.734846711099164,1e4), ...
4.                   @(t)mat4init(t,g,v,theta),k);
5. opts_jac = bvpset('FJacobian',@(t,u,k)jac(t,u,k,n), ...
6.              'BCJacobian',@(u0,ut,k)bc_ab(u0,ut,k), ...
7.              'RelTol',0.1,'AbsTol',0.1,'Stats','on');
8. opts = bvpset('RelTol',0.1,'AbsTol',0.1,'Stats','on');
9. tic
10. sol = bvp4c(  @(t,u,k)mat4ode(t,u,k,n,g), ...
11.             @(u0,ut,k)mat4bc(u0,ut,k,v,theta), solinit,opts_jac);
12. toc
13. % 输出验证
14. fprintf('\n阻力系数k=：%7.3f.\n',sol.parameters)
15. t=linspace(0,0.7348,1e3)';
16. ut = deval(sol,t)';
17. x=ut(:,1);y=ut(:,2);
18. plot(x,y)
19. % 微分方程
20. function du = mat4ode(~,u,k,n,g)
21. vx = u(3); vy = u(4);
22. v  = sqrt(vx^2+vy^2);
23. dx  = vx;
24. dy  = vy;
25. dvx =   -k*v^(n-1)*vx;
26. dvy = -g-k*v^(n-1)*vy;
27. du  = [dx dy dvx dvy]';
28. end
29. % FJacobian
30. function [J,du_k]= jac(~,u,k,n)
31. vx = u(3); vy = u(4);
32. v  = sqrt(vx^2+vy^2);
33. dvx_vx = -k*v^(n-1)-k*(n-1)*v^(n-3)*vx*vx;
34. dvx_vy = -k*(n-1)*v^(n-3)*vx*vy;
35. dvy_vx = -k*(n-1)*v^(n-3)*vx*vy;
36. dvy_vy = -k*v^(n-1)-k*(n-1)*v^(n-3)*vy*vy;
37. J = [
38.     0 0 1 0;
39.     0 0 0 1;
40.     0 0 dvx_vx dvx_vy;
41.     0 0 dvy_vx dvy_vy];
42. du_k=[0 0 -v^(n-1)*vx -v^(n-1)*vy]';
43. end
44. % 边界条件
45. function bc = mat4bc(u0,ut,~,v,theta)
46. x0=u0(1);y0=u0(2);vx0=u0(3);vy0=u0(4);
47. xt=ut(1);yt=ut(2);vxt=ut(3);vyt=ut(4);
48. ua=[x0 y0 vx0-v*cos(theta) vy0-v*sin(theta)]';
49. bc = [ua(1:4);yt];
50. end
51. % BCJacobian
52. function [a,b,k]=bc_ab(~,~,~)
53. a=[ 1 0 0 0 0
54.     0 1 0 0 0
55.     0 0 1 0 0
56.     0 0 0 1 0
57.     0 0 0 0 0 ];
58. b=[ 0 0 0 0 0
59.     0 0 0 0 0
60.     0 0 0 0 0
61.     0 0 0 0 0
62.     0 1 0 0 0 ];
63. k=[ 0 0 0 0 0 ]';
64. end
65. % 初值估计
66. function u_t = mat4init(t,g,v,theta)
67. vxt=v*cos(theta);vyt=v*sin(theta);
68. u_t=[vxt*t vyt*t-0.5*g*t^2 vxt vyt]';
69. end