Advertisement
Y-BH

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
Advertisement
RAW Paste Data Copied
Advertisement