用Matlab改進一次二階矩法程序

1、用Matlab編的計算結構可靠指標的改進一次二階矩法程序（驗算點法） 題目：編制改進一次二階矩法計算可靠指標的程序，并給出算例，要求提供源程序選取的算例為：z=g(x,y)=x*y-1140，其中x,y服從正態(tài)分布，x=38,Vx=0.1, y=38,Vy=0.05本程序采用Matlab編寫。選取1=3.0，2=2.5計算結果為：可靠指標4.2672，最終驗算點為：(22.8430 , 49.9060)，在驗算點處功能函數(shù)值為：1.2354e-004%保存為strRlbt.m，在Matlab命令窗口中輸入strRlbt執(zhí)行即可N = 2;%變量個數(shù)miu = 38 54;%均值v = 0.1

2、0.05;%變異系數(shù)sgma = miu .* v;%方差syms x yg = sym(x * y - 1140);%功能函數(shù)jacg = jacobian( g ,x;y);%計算雅可比矩陣initvalue = miu;v;sgma;%用作函數(shù)參數(shù)%選取beta，定義x0=miubeta1 = 3.0;xopt0 = 38 54;alpha0 = zeros(1,2); alpha1 , xopt1 , result = calforbeta( initvalue , beta1 , alpha0 , xopt0 , jacg , g );if result = 1 disp(第一次假定

3、的可靠指標); returnend %再次假定betabeta2 = 2.5;xopt0 = miu - beta2 * alpha1 .* sgma;gvalue = jacgfunc(jacg,xopt1);alpha0 = (sgma .* gvalue) / sqrt(sum(sgma .* gvalue).2); alpha2 , xopt2 , result = calforbeta( initvalue , beta2 , alpha0 , xopt0 , jacg , g );if result = 1 disp(第二次假定的可靠指標); returnend%beta迭代求解g

4、1 = gfunc(g,xopt1);g2 = gfunc(g,xopt2);eps = 0.1; %精度while abs(g2) eps temp = beta2; beta2 = beta2 - (beta2 - beta1)/(g2 - g1) * g2; beta1 = temp; alpha2 , xopt2 , result = calforbeta( initvalue , beta2 , alpha1 , xopt1 , jacg , g ); temp = g2; g2 = gfunc(g,xopt2); g1 = temp; if result = 1 break end

5、enddisp(可靠指標為：);disp(beta2);disp(最終驗算點為：);disp(xopt2);disp(在驗算點處功能函數(shù)值為：);disp(g2);function g_out = gfunc( g , x_in )%功能函數(shù)值計算x = x_in(1);y = x_in(2);g_out = eval(g);%函數(shù)值%將以上內(nèi)容保存為gfunc.mfunction g_out = jacgfunc( jacg , x_in )%功能函數(shù)偏導數(shù)計算，即雅可比矩陣計算x = x_in(1);y = x_in(2);for i = 1:2 g_out(i) = eval(jacg(

6、i);%1為對x的導數(shù),2為對y的導數(shù)end%將以上內(nèi)容保存為jacgfunc.mfunction alpha1 , xopt1 ,result = calforbeta( initvalue , beta0 , alpha0 , xopt0 , jacg , g)%對選取的beta進行計算result = 0;N = length(xopt0);alpha = alpha0;xopt = xopt0;%initvalue為初始值miu = initvalue(1,:);%第一行為均值v = initvalue(2,:);%第二行為變異系數(shù)sgma = initvalue(3,:);%第三行為

7、方差eps = 0.1;while 1 %功能函數(shù)達到精度則退出循環(huán),result=1表示計算出可靠指標 if abs(gfunc(g,xopt0) eps alpha1 = alpha0; xopt1 = xopt0; result = 1; break; end %計算alpha和新的驗算點xopt gvalue = jacgfunc(jacg,xopt); sgmaz = sqrt(sum(sgma .* gvalue).2); alpha0 = sgma .* gvalue / sgmaz; xopt0 = miu - beta0 * alpha0 .* sgma; sum1 = sum(alpha - alpha0).2); sum2 = sum(xopt - xopt0).2); alpha =