一次二阶矩程序

civilmanhuijun

请问哪位有一次二阶矩的源程序啊？谢谢了，小弟进来急需，再次感谢了
vb，fortran或matlab都可以。

风花雪月

一个例子:

1. 题目：编制改进一次二阶矩法计算可靠指标的程序，并给出算例，要求提供源程序，算法语言不限。
   
2. 
   
3. 选取的算例为：z=g(x,y)=x*y-1140，其中x,y服从正态分布，μx=38,Vx=0.1, μy=38,Vy=0.05
   
4. 
   
5. 本程序采用Matlab编写。
   
6. 
   
7. 选取β1=3.0，β2=2.5
   
8. 
   
9. 计算结果为：可靠指标β＝4.2672，最终验算点为：(22.8430 , 49.9060)，在验算点处功能函数值为：1.2354e-004
   
10. 
    
11. 
    
12. 
    
13. %保存为strRlbt.m，在Matlab命令窗口中输入strRlbt执行即可
    
14. 
    
15. N = 2;%变量个数
    
16. 
    
17. miu = [38 54];%均值
    
18. 
    
19. v = [0.1 0.05];%变异系数
    
20. 
    
21. sgma = miu .* v;%方差
    
22. 
    
23. syms x y
    
24. 
    
25. g = sym('x * y - 1140');%功能函数
    
26. 
    
27. jacg = jacobian( g ,[x;y]);%计算雅可比矩阵
    
28. 
    
29. initvalue = [miu;v;sgma];%用作函数参数
    
30. 
    
31. %选取beta，定义x0=miu
    
32. 
    
33. beta1 = 3.0;
    
34. 
    
35. xopt0 = [38 54];
    
36. 
    
37. alpha0 = zeros(1,2);
    
38. 
    
39. [ alpha1 , xopt1 , result ] = calforbeta( initvalue , beta1 , alpha0 , xopt0 , jacg , g );
    
40. 
    
41. if result == 1
    
42. 
    
43.    disp('第一次假定的饧次煽恐副辏');
    
44. 
    
45.    return
    
46. 
    
47. end
    
48. 
    
49. 
    
50. 
    
51. %再次假定beta
    
52. 
    
53. beta2 = 2.5;
    
54. 
    
55. xopt0 = miu - beta2 * alpha1 .* sgma;
    
56. 
    
57. gvalue = jacgfunc(jacg,xopt1);
    
58. 
    
59. alpha0 = (sgma .* gvalue) / sqrt(sum((sgma .* gvalue).^2));
    
60. 
    
61. [ alpha2 , xopt2 , result ] = calforbeta( initvalue , beta2 , alpha0 ,  xopt0 , jacg , g );
    
62. 
    
63. if result == 1
    
64. 
    
65.    disp('第二次假定的饧次煽恐副辏');
    
66. 
    
67.    return
    
68. 
    
69. end
    
70. 
    
71. 
    
72. 
    
73. %beta迭代求解
    
74. 
    
75. g1 = gfunc(g,xopt1);
    
76. 
    
77. g2 = gfunc(g,xopt2);
    
78. 
    
79. eps = 0.1; %精度
    
80. 
    
81. while abs(g2) > eps
    
82. 
    
83.    temp = beta2;
    
84. 
    
85.    beta2 = beta2 - (beta2 - beta1)/(g2 - g1) * g2;
    
86. 
    
87.    beta1 = temp;
    
88. 
    
89.    [ alpha2 , xopt2 , result ] = calforbeta( initvalue , beta2 , alpha1 ,  xopt1 , jacg , g );
    
90. 
    
91.    temp = g2;
    
92. 
    
93.    g2 = gfunc(g,xopt2);
    
94. 
    
95.    g1 = temp;
    
96. 
    
97.    if result == 1
    
98. 
    
99.       break
    
100. 
     
101.    end
     
102. 
     
103. end
     
104. 
     
105. disp('可靠指标为：');
     
106. 
     
107. disp(beta2);
     
108. 
     
109. disp('最终验算点为：');
     
110. 
     
111. disp(xopt2);
     
112. 
     
113. disp('在验算点处功能函数值为：');
     
114. 
     
115. disp(g2);
     
116. 
     
117. 
     
118. 
     
119. function g_out = gfunc( g , x_in )
     
120. 
     
121. %功能函数值计算
     
122. 
     
123. x = x_in(1);
     
124. 
     
125. y = x_in(2);
     
126. 
     
127. g_out = eval(g);%函数值
     
128. 
     
129. %将以上内容保存为gfunc.m
     
130. 
     
131. 
     
132. 
     
133. function g_out = jacgfunc( jacg , x_in )
     
134. 
     
135. %功能函数偏导数计算，即雅可比矩阵计算
     
136. 
     
137. x = x_in(1);
     
138. 
     
139. y = x_in(2);
     
140. 
     
141. for i = 1:2
     
142. 
     
143.    g_out(i) = eval(jacg(i));%1为对x的导数,2为对y的导数
     
144. 
     
145. end
     
146. 
     
147. %将以上内容保存为jacgfunc.m
     
148. 
     
149. 
     
150. 
     
151. function [ alpha1 , xopt1 ,result ] = calforbeta( initvalue , beta0 , alpha0 , xopt0 , jacg , g)
     
152. 
     
153. %对选取的beta进行计算
     
154. 
     
155. result = 0;
     
156. 
     
157. N = length(xopt0);
     
158. 
     
159. alpha = alpha0;
     
160. 
     
161. xopt = xopt0;
     
162. 
     
163. %initvalue为初始值
     
164. 
     
165. miu = initvalue(1,:);%第一行为均值
     
166. 
     
167. v = initvalue(2,:);%第二行为变异系数
     
168. 
     
169. sgma = initvalue(3,:);%第三行为方差
     
170. 
     
171. eps = 0.1;
     
172. 
     
173. while 1
     
174. 
     
175.    %功能函数达到精度则退出循环,result=1表示计算出可靠指标
     
176. 
     
177.    if abs(gfunc(g,xopt0)) < eps
     
178. 
     
179.       alpha1 = alpha0;
     
180. 
     
181.       xopt1 = xopt0;
     
182. 
     
183.       result = 1;
     
184. 
     
185.       break;
     
186. 
     
187.    end
     
188. 
     
189.    %计算alpha和新的验算点xopt
     
190. 
     
191.    gvalue = jacgfunc(jacg,xopt);
     
192. 
     
193.    sgmaz = sqrt(sum((sgma .* gvalue).^2));
     
194. 
     
195.    alpha0 = sgma .* gvalue / sgmaz;
     
196. 
     
197.    xopt0 = miu - beta0 * alpha0 .* sgma;
     
198. 
     
199.    sum1 = sum((alpha - alpha0).^2);
     
200. 
     
201.    sum2 = sum((xopt - xopt0).^2);
     
202. 
     
203.    alpha = alpha0;
     
204. 
     
205.    xopt = xopt0;
     
206. 
     
207.    %醝和验算点xi达到精度则退出循环
     
208. 
     
209.    if sum1 < 0.001 | sum2 < 0.001
     
210. 
     
211.       alpha1 = alpha0;
     
212. 
     
213.       xopt1 = xopt0;
     
214. 
     
215.       break;
     
216. 
     
217.    end
     
218. 
     
219. end
     
220. 
     
221. %将以上内容保存为calforbeta.m

复制代码

来自http://littlebees.blog.hexun.com/14143949_d.html

civilmanhuijun

非常感谢风花雪月的热情回复，谢谢
收下了:lol

Tony1027

第一次假定的饧次煽恐副辏？什么意思？。。。