分享一个打靶法求解两点边值问题的Matlab函数

Catsayer

1. function varargout = shooting_two_point_boundary(varargin)
   
2. % ==========================================================
   
3. % 函数名：shooting_two_point_boundary.m
   
4. % 基于打靶法计算两点边值问题，仅针对二阶微分方程
   
5. % author: xianfa110.
   
6. % blog： <a href="http://blog.sina.com.cn/xianfa110" target="_blank">http://blog.sina.com.cn/xianfa110</a>
   
7. % 函数形式：
   
8. % [result,err,z0] = shooting_two_point_boundary(@fun,[y_0,y_end],[x_0,x_1],h);
   
9. % 输入：
   
10. % fun = 函数名;
    
11. % y_0 = 函数初值;
    
12. % y_end = 函数终值;
    
13. % x_0 = 自变量初值;
    
14. % x_end = 自变量终值;
    
15. % h = 积分步长;
    
16. % 输出：
    
17. % result = [x,y];
    
18. % err = 误差;
    
19. % z0 = y'初值;
    
20. % ===========================================================
    
21. % 函数fun：4y''+yy' = 2x^3 +16 ; 2<= x <=3
    
22. % 写法：
    
23. % function f = fun(y,x)
    
24. % dy = y(2);
    
25. % dz = (2*x^3+16-y(1)*y(2))/4;
    
26. % f = [dy,dz];
    
27. % ===========================================================
    
28. % 注意:y(1) = y，y(2) = y'。
    
29. % ===========================================================
    
30. F = varargin{1};
    
31. y_0 = varargin{2}(1);
    
32. y_end = varargin{2}(2);
    
33. x_0 = varargin{3}(1);
    
34. x_1 = varargin{3}(2);
    
35. ts = varargin{4};
    
36. t0 = x_0-0.5;
    
37. flg = 0;
    
38. kesi = 1e-6;
    
39. y0 = rkkt(F,[y_0,t0],x_0,x_1,ts);
    
40. n = length(y0(:,1));
    
41. if abs(y0(n,1)-y_end)<=kesi
    
42. flg=1;
    
43. else
    
44. t1=t0+1;
    
45. y1=rkkt(F,[y_0,t1],x_0,x_1,ts);
    
46. if abs(y1(n,1)-y_end)<=kesi
    
47. flg=1;
    
48. end
    
49. end
    
50. if flg ~= 1
    
51. while abs(y1(n,1)-y_end) > kesi
    
52. % ==========插值法求解非线性方程=============== %
    
53. t2 = t1-(y1(n,1)-y_end)*(t1-t0)/(y1(n,1)-y0(n,1));
    
54. y2 = rkkt(F,[y_0,t2],x_0,x_1,ts);
    
55. t0=t1;
    
56. t1=t2;
    
57. y0=y1;
    
58. y1=y2;
    
59. end
    
60. end
    
61. x = x_0:ts:x_1;
    
62. out = [x',y1(:,1)];
    
63. varargout{1} = out;
    
64. varargout{2} = abs(y1(n,1)-y_end);
    
65. varargout{3} = t1;

复制代码

参考文献：数值计算原理
相关文章：
自编基于龙格库塔法的Matlab数值积分函数 (本文所用的积分函数)
Matlab非线性方程求解器fsolve总结(含实例)

转自：http://blog.sina.com.cn/s/blog_408540af0100b7mi.html