MATLAB解决一次二阶矩问题

1.功能函数为：

2.变量分布为

变量X1、X2服从的分布为N(10,5)的正态分布

允许的设计精度为：

 3.源代码段

%%假定随机变量X1和随机变量X2服从正态分布，验算点法计算可靠度指标B，允许误差设置在1e-3
%流程如下：
%1.假设初始设计点，一般设置在均值点处
%2.计算出失效概率
%3.计算出灵敏度系数（cosx1和cosx2）
%4.计算出新的设计点坐标值
%5.判断是否满足精度，若不满足继续迭代
%6.验算生成的最终设计点是否符合限制条件
%%初始化参数
clc;clear
syms x1 x2; 
Z=x1^3+x2^3-18          		%这里设置的是功能函数
dzx1=diff (Z,'x1')     			%求功能函数Z对x1的偏导数
dzx2=diff (Z,'x2')      		%求功能函数Z对x2的偏导数
Z1=[dzx1,dzx2]         			%Z1为梯度矩阵（数据类型为sym）     
x=[10,10]           		    %变量1和变量2服从的分布的均值
sigmax=[5,5]        			%设置的是变量服从分布的标准差

origin_x=[10,10]
%% 迭代部分
x0=repmat(eps,length(x),1);        				    %注：eps为计算机识别最小数，进入循环
i=0;
while norm(x-x0)>1e-3              				    %用while来控制解的精度
    x0=x;                                   		%将x赋给x0，方便后续进行变量对比
    i=1+i;
    x1=x0(1)                						%定义全局变量x1，初始化点为均值点
    x2=x0(2)                						%定义全局变量x2，初始化点为均值点
    grad=double(subs(Z1))   					%将梯度矩阵转换为double类型 
    miuzl=x(1)^3+x(2)^3-18+grad*(origin_x-x)' 	%计算功能函数均值,这里为化简式。
    sigmaxzl=norm(sigmax.*grad)   				%计算功能函数的标准差         
    B=miuzl/sigmaxzl                      		%计算出B（失效概率）
    cosx1=-(sigmax(1)*grad(1))/(sigmaxzl)     	%计算出灵敏度系数1 
    cosx2=-(sigmax(2)*grad(2))/(sigmaxzl)     	%计算出灵敏度系数2 
    newx1=origin_x(1)+sigmax(1)*B*cosx1
    newx2=origin_x(2)+sigmax(2)*B*cosx2
    x=[newx1,newx2]                      		%差值是否满足误差？若不满足接着迭代
    B_data(i)=miuzl/sigmaxzl
    X1_data(i,:)=x                              
end 
%% 验算部分
final_Z=x1^3+x2^3-18;                           %对最后收敛的设计点进行验算
if final_Z<1e-3
    fprintf('可靠度指标β=%f，',B)
    fprintf('设计点横坐标为%f，',x(1))
    fprintf('设计点纵坐标为%f，且设计点满足限制条件。',x(2))
end