echo off
clear all;
clc;
%定义全局变量；
global wr wa wb T_mr v_t vi;
%v=zeros(1,5);
rho=1.29;         %空气密度 单位为kg/m^3 
R=1.5;            %旋翼半径  单位为m
Omega=870*2*pi/60;% 主旋翼角速度  单位为rad/s
b1=1/180;         %纵向周期变距角  单位为rad
is=5/180;         %桨盘倾角  单位为rad
a=6.0;            %升力曲线斜率 单位为1/rad 
b=2;              %桨叶数
c=0.354*0.3048;   % 平均桨叶弦长[m] 
theta_coll=10/180;%总距角  单位为rad
theta_twist=-5/180;%负扭  单位为rad
ua=10;             %x轴线速度  单位为m/s
va=10;             %y轴线速度  单位为m/s
wa=10;             %z轴线速度  单位为m/s
A=pi*R^2;          %桨盘面积   单位为m^2
a1=1/180;         %横向周期变距  单位为rad
%for a1=-10/180:1/180:10/180;    %横向周期变距角一组变化值  单位为rad
    %wr=zeros(1,10);
    %wb=zeros(1,10);
    wr=wa+(a1+is)*ua-b1*va;
    wb=wr+2/3*R*(theta_coll+3/4*theta_twist);
    for i=1:5;
    %i=1;
    v(i)=0;
    T_mr=(wb-v(i))*Omega*R^2*rho*a*b*c/4;
    %v_t=zeros(1,10);
    %vi=zeros(1,10);
    v_t=(ua^2+va^2+wa^2*(wr-2*v(i)))^(1/2);
    v1_1=((1/4*v_t^4+(T_mr/(2*rho*pi*R*R))^2)^(1/2)-1/2*v_t^2)^(1/2);
    if abs(v1_1-v(i))<eps
        v(i)=v1_1;
    break
    else i=i+1;
        v(i)=v1_1;
    %G(vi)=1/4*v_t^4+T_mr^2/4*rho^2*A^2;  %任意假设函数，为了方便计算
    %H(vi)=1/2*v_t^2;                     %任意假设函数，为了方便计算
    %f(vi)=vi^2-G(vi)^1/2-H(vi)
    %vi=fsolve(@f(vi) 0);
    end
  end
%end
    hold on;
 plot(a1,v(i));
 xlabel('纵向周期变距a1（单位rad）');
 ylabel('诱导速度vi（单位m/s）');
 title('诱导速度随纵向周期变距变化图');
 text(5/180,v(i));
 grid on
程序大致就是要执行压缩的图中的一些方程的迭代，但是怎么写都不对，希望朋友们能帮我看看错在哪了，谢谢了，论文写不出来，程序也编不出来，急疯了！