精华内容
下载资源
问答
  • 我们要求的就是theta角和时间t之间的关系曲线,这是一道典型的二阶常微分方程的求解,用四阶龙格库塔方程可以求解。 2 古典龙格库塔算法公式: 高等数值计算课本(清华大学出版社,186页) 但是古典龙格库塔方法...

    前几天接到师姐分派的任务,让我求解一艘船模的横摇角的时间历程曲线,为后期的减摇控制准备。

    1 首先冷静分析一下,原方程如下:

    在这里插入图片描述
    我们要求解的就是theta角和时间t之间的关系曲线,这是一道典型的二阶常微分方程的求解,用四阶龙格库塔方程可以求解。

    2 古典龙格库塔算法公式:

    节选自高等数值计算
    高等数值计算课本(清华大学出版社,186页)

    但是古典龙格库塔方法解决的是一阶常微分方程,也就是类似于这样的方程组

    在这里插入图片描述
    但是我们现在面对的是
    在这里插入图片描述
    查了一下书,没有教二阶常微分该怎么解,于是自己网上找答案:

    3 网上的一道例题提供了解决办法

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    4 现在可以开始编程了,我没有用Matlab自带的Ode45函数,我是自己编写的四阶Rungekutta方程,代码如下:

    %横摇角时间历程曲线计算主程序
    global a b c d h t y z D h_
    a = 18.4;
    b = 0.3;
    c = 214.88;
    d = 0;
    h = 0.1;
    y = 0;
    z = 0;
    D = 1.58*10^3;
    h_ = 0.136;
    t = 0;
    for i = 1:1000
       
       d = -D*h_*0.26;
        d1=sin(t);
        
        [Y,Z] = RungeKutta_2(y,z);
        t = t+h;
        z = Z;
        y = Y;
        y_(i) = y
        t_(i)=t;
    end
    max(y)
    plot(t_(1:1000),y_(1:1000));
    xlabel('时间t')
    ylabel('横摇角theta')
    title('横摇角时历曲线')
    

    该程序引用了Rungekutta_2.m程序,代码如下:

    function [Y,Z]=RungeKutta_2(y,z)
    global a b c d h y z D h_ t 
    
    k11=z;
    k21=(d*sin(t)-c*y-b*z)/a;
    k12=z+(h/2)*k21;
    k22=(d*sin(t+0.5*h)-c*(y+(h/2)*k11)-b*(z+(h/2)*k21))/a;
    k13=z+(h/2)*k22;
    k23=(d*sin(t+0.5*h)-c*(y+(h/2)*k12)-b*(z+(h/2)*k22))/a;
    k14=z+h*k23;
    k24=(d*sin(t+h)-c*(y+h*k13)-b*(z+h*k23))/a;
    
    Z=z+(h/6)*(k11+2*k12+2*k13+k14);    
    Y=y+(h/6)*(k21+2*k22+2*k23+k24);
    

    用的时候把这两个放在一个文件夹就行了。

    5 运行结果

    在这里插入图片描述
    符合规则波下船舶的横摇运动特征。完事了。

    一起学习,一起进步!!

    展开全文
  • 方程求解程序清单 a=-1,b=2,c=-1; w=1; m=2; n=1; h = 0.02; t=0:h:30; s1=dsolve('a*D2y+b*Dy+c*y=sin(w*t)','y(0)=m,Dy(0)=n','t'); s1_n = eval(s1); hold on plot(t,s1_n,'ko'); EulerOED(a,b,c,w,m,n,h); hold ...

    方程求解程序清单 a=-1,b=2,c=-1; w=1; m=2; n=1; h = 0.02; t=0:h:30; s1=dsolve('a*D2y+b*Dy+c*y=sin(w*t)','y(0)=m,Dy(0)=n','t'); s1_n = eval(s1); hold on plot(t,s1_n,'ko'); EulerOED(a,b,c,w,m,n,h); hold off function EulerOED(a,b,c,w,x0,x1,h) A = [x0;x1]; t=0:h:30; for i = 1:1:length(t)-1 A(:,i+1) = [1,h;(-(c/a)*h),(1-(b/a)*h)]*A(:,i) + [0;(h/a)]*sin(w*t(i)); end plot(t,A(1,:),'r*'); 对于二阶全微分方程a*y''(t)+b*y'(t)+c=sin(wt) ,不同的a,b,c,w取值和初始条件会求出不同的解,通解又是由齐次解和特解组成。其中,齐次解由特征方程决定,而特解的决定因素则比较复杂。 讨论思路 (1)通解随初始条件变化情况 (2)通解随a,b,c变化情况 b^2-4ac>0(两个不同的实根) b^2-4ac=0(两个相同的重根) b^2-4ac<0(两个不同的复数根) 1).b>0 2).b=0 3).b<0 (3)通解随w变化情况 b^2-4ac=0情况 b^2-4ac<0情况 (3)通解随w变化的规律 W属于(0,1)时,随w的增大在齐次解的旁边波动 w属于(1,+),随w的增大逐渐趋近于齐次解。 Matlab解二阶常微分方程 方程:a*y''(t)+b*y'(t)+c=sin(wt) 求解:1.解析解 2.数值解(欧拉方法) 目的:1.比较两种求解方式的拟合情况 2.通解随w变化的规律 (1)通解随初始条件变化情况 Ex1: a=2,b=3,c=1,y(0)=0;y'(0)=0,w=1 Ex2: a=2,b=3,c=1,y(0)=2;y'(0)=0,w=1 Ex3: a=2,b=3,c=1,y(0)=2;y'(0)=4,w=1 (2)通解随a,b,c变化情况 Ex1: a=2,b=3,c=1,y(0)=0;y'(0)=0,w=1 Ex4: a=-2,b=3,c=1,y(0)=0;y'(0)=0,w=1 Ex5: a=2,b=-3,c=1,y(0)=0;y'(0)=0,w=1 Ex6: a=2,b=3,c=-1,y(0)=2y'(0)=1,w=1 ? EX: a=2 ,b=2*sqrt(2) ,c=1,y(0)=0;y'(0)=0,w=1 (3).b^2-4ac<0 EX:a=4,b=-1,c=2,y(0)=0;y'(0)=0,w=1 EX:a=4,b=1,c=2,y(0)=3,y'(0)=0,w=1 EX:a=4,b=0,c=1,y(0)=2;y'(0)=0,w=1

    展开全文
  • 今年年初的时候给师姐做了DDPG算法的船舶减横摇控制算法,师姐还有想法要让我把纵摇-埀荡两个自由度的减摇也做出来,这个任务归我了。实际上不管是多少个自由...这个方程组是典型的二阶常微分方程组,一堆的水动力参数

    今年年初的时候给师姐做了DDPG算法的船舶减横摇控制算法,师姐还有想法要让我把纵摇-埀荡两个自由度的减摇也做出来,这个任务归我了。实际上不管是多少个自由度的减摇,其实都需要解运动方程,当初做单自由度横摇的时候就是用四阶龙格库塔求解二阶微分方程,具体代码见:

    https://blog.csdn.net/Mezikov/article/details/103898387

    0 写在前面

    代码内容在第二节,第一节是运动方程各种参数求解。

    1 冷静分析一下方程组

    在这里插入图片描述
    这个方程组是典型的二阶常微分方程组,一堆的水动力参数,但都是可以求出来的已知量(求解方法看下图,Z和M分别是方程1等号右边项和方程2等号右边项),未知量是埀荡量ZG和纵摇角theta及其一阶二阶导。
    在这里插入图片描述
    在这里插入图片描述

    2 现在可以开始编程了,我没有用Matlab自带的Ode45函数,我是自己编写的四阶Rungekutta方程,代码如下:

    首先是纵摇运动方程 pitch equation

    function w2 = pitch_equation_test3(y1,z1,y2,z2)
    % global B d delta V v g L rho A_w C_w C_p C_b W_e x_b GM_l b_zz fi
    fi = pi/3;       % 波浪遭遇角
    B = 1;           % width of boat
    d = 0.341;         % draft
    delta = 1.76*10^3;  % displacement
    V = 1.76;           % displacement volume
    g = 9.8;
    L = 6;             % length
    lamda = 6;         % 波长
    rho = 1*10^3;
    A_w = 5.448;        % area of water plane
    A_m = 0.346;        % '中横剖面面积'
    C_w = A_w / (L * B);  % coefficient of water plane
    C_p = V/(L*A_m);
    C_b = delta/(1000*L*B*d);
    v = 4.44;              %船速
    T = 0.8*sqrt(lamda);
    W_0 = 2*pi/T  ;      % 波浪圆频率
    W_e = W_0*(1+v*W_0*cos(fi)/g) ;   % 波浪遭遇频率'%%%%%%%%%%%%%%%%%
    x_b = -0.23 ;     % 'center of float' 这个很可能是错的,因为我自己变得,按照大尺度船的计算是0.23,但是根据参考文献,这里应该是负的
    GM_l = (L^2*(5.55*C_w+1)^3/3450)/(C_b*d);
    M = 1.6*10^3;      % weight
    H_0 = B/(2*d);
    a_zz = 0.8*H_0*C_w*M;
    b_zz = (5.4*(C_w/C_p)*(H_0)^(0.5)-4.7)*delta/((g*L)^0.5);
    c_zz = rho*g*A_w;
    a_z0 = (v/W_e^2)*b_zz-a_zz*x_b;
    b_z0 = -b_zz*x_b+v*a_zz;     % 0其实就是theta
    c_z0 = -rho*g*A_w*x_b;
    I_00 = 0.07*C_w*M*L^2;
    a_00 = 0.83*H_0*C_p^2*(0.25*L)^2*(delta/g)^0.5;
    b_00 = 0.08*H_0*delta*L^2/(g*L)^0.5;
    c_00 = rho*g*(V*GM_l);
    a_0z = -(a_zz*x_b);
    b_0z = -(b_zz*x_b);
    c_0z = -rho*g*A_w*x_b-v*b_zz;
    fi = pi/3;
    a = 0.5*0.17*lamda^(3/4);
    k = 2*pi/lamda;
    n1 = (I_00+a_00)-(a_0z*a_z0)/(M+a_zz);
    n2 = rho*g*a*k*exp(-k*d/2)*B*d*sin(k*B*L*sin(fi)/2)/(k*B*L*sin(fi)/2);
    n3 = (I_00+a_00)*(M+a_zz)-a_0z*a_z0;
      
    w2 =  (b_00*(M+a_zz)-a_0z*b_z0)/n3*z2+...
          (c_00*(M+a_zz)-a_0z*c_z0)/n3*y2+...
          (b_0z*(M+a_zz)-a_0z*b_zz)/n3*z1+...
          (c_0z*(M+a_zz)-a_0z*c_zz)/n3*y1;
                 
    end
    

    然后是埀荡运动方程 heave equation

    function w1 = heave_equation_test3(y1,z1,y2,z2)
    % global B d delta V v g L rho A_w C_w C_p C_b W_e x_b GM_l b_zz fi
    fi = pi/3;       % 波浪遭遇角2*pi/3
    B = 1;           % width of boat
    d = 0.341;         % draft
    delta = 1.76*10^3;  % displacement
    V = 1.76;           % displacement volume
    g = 9.8;
    L = 6;             % length
    lamda = 6;         % 波长
    rho = 1*10^3;
    A_w = 5.448;        % area of water plane
    A_m = 0.346;        % '中横剖面面积'
    C_w = A_w / (L * B);  % coefficient of water plane
    C_p = V/(L*A_m);
    C_b = delta/(1000*L*B*d);
    v = 4.44;              % 这个是我瞎编的
    T = 0.8*sqrt(lamda);
    W_0 = 2*pi/T ;      % 波浪圆频率
    W_e = W_0*(1+v*W_0*cos(fi)/g) ;     % 波浪遭遇频率'%%%%%%%%%%%%%%%%%
    x_b = -0.23   ;     % 'center of float' 这个很可能是错的,因为我自己变得,按照大尺度船的计算是0.23,但是根据参考文献,这里应该是负的
    GM_l = (L^2*(5.55*C_w+1)^3/3450)/(C_b*d);
    M = 1.76*10^3;      % weight
    H_0 = B/(2*d);
    a_zz = 0.8*H_0*C_w*M;
    b_zz = (5.4*(C_w/C_p)*(H_0^0.5)-4.7)*delta/((g*L)^0.5);
    c_zz = rho*g*A_w;
    a_z0 = (v/W_e^2)*b_zz-a_zz*x_b;
    b_z0 = -b_zz*x_b+v*a_zz;     % 0其实就是theta
    c_z0 = -rho*g*A_w*x_b;
    I_00 = 0.07*C_w*M*L^2;
    a_00 = 0.83*H_0*C_p^2*(0.25*L)^2*(delta/g)^0.5;
    b_00 = 0.08*H_0*delta*L^2/(g*L)^0.5;
    c_00 = rho*g*(V*GM_l);
    a_0z = -(a_zz*x_b);
    b_0z = -(b_zz*x_b);
    c_0z = -rho*g*A_w*x_b-v*b_zz;
    fi =pi/3;
    a = 0.5*0.17*lamda^(3/4);
    k = 2*pi/lamda;
    n1 = (I_00+a_00)-(a_0z*a_z0)/(M+a_zz);
    n2 = rho*g*a*k*exp(-k*d/2)*B*d*sin(k*B*L*sin(fi)/2)/(k*B*L*sin(fi)/2);
    n3 = (I_00+a_00)*(M+a_zz)-a_0z*a_z0;
    
    w1 =  (b_zz*(I_00+a_00)-a_z0*b_0z)/n3*z1+...
          (c_zz*(I_00+a_00)-a_z0*c_0z)/n3*y1+...
          (b_z0*(I_00+a_00)-a_z0*b_00)/n3*z2+...
          (c_z0*(I_00+a_00)-a_z0*c_00)/n3*y2;
    
    end
    

    然后是我的龙格库塔算法代码Rungukutta_3

    clear
    clc
    
    y1(1)=0.1;
    z1(1)=0.1;
    y2(1)=0.1;
    z2(1)=0.04;
    h=0.1;
    t=0:h:500;
    for i=1:length(t)-1 
        
        k1 = z1(i); 
        k1_= z2(i);  
        L1 = cal_F(t(i))-heave_equation_test3(y1(i),z1(i),y2(i),z2(i));
        L1_= cal_M(t(i))-pitch_equation_test3(y1(i),z1(i),y2(i),z2(i));
        k2 = z1(i)+0.5*h*L1;
        k2_= z2(i)+0.5*h*L1_;
        L2 = (cal_F(t(i))+cal_F(t(i)+h))/2-heave_equation_test3(y1(i)+0.5*h*k1,z1(i)+0.5*h*L1,y2(i)+0.5*h*k1_,z2(i)+0.5*h*L1_);
        L2_= (cal_M(t(i))+cal_M(t(i)+h))/2-pitch_equation_test3(y1(i)+0.5*h*k1,z1(i)+0.5*h*L1,y2(i)+0.5*h*k1_,z2(i)+0.5*h*L1_);
        k3 = z1(i)+0.5*h*L2;
        k3_= z2(i)+0.5*h*L2_;
        L3 = (cal_F(t(i))+cal_F(t(i)+h))/2-heave_equation_test3(y1(i)+0.5*h*k2,z1(i)+0.5*h*L2,y2(i)+0.5*h*k2_,z2(i)+0.5*h*L2_);
        L3_= (cal_M(t(i))+cal_M(t(i)+h))/2-pitch_equation_test3(y1(i)+0.5*h*k2,z1(i)+0.5*h*L2,y2(i)+0.5*h*k2_,z2(i)+0.5*h*L2_);
        k4 = z1(i)+h*L3;
        k4_= z2(i)+h*L3_;
        L4 = (cal_F(t(i))+cal_F(t(i)+h))/2-heave_equation_test3(y1(i)+h*k3,z1(i)+h*L3,y2(i)+h*k3_,z2(i)+h*L3_);
        L4_= (cal_M(t(i))+cal_M(t(i)+h))/2-pitch_equation_test3(y1(i)+h*k3,z1(i)+h*L3,y2(i)+h*k3_,z2(i)+h*L3_);
        y1(i+1)=(y1(i)+(h/6)*(k1+2*k2+2*k3+k4));
        z1(i+1)=(z1(i)+(h/6)*(L1+2*L2+2*L3+L4));
        y2(i+1)=(y2(i)+(h/6)*(k1_+2*k2_+2*k3_+k4_));
        z2(i+1)=(z2(i)+(h/6)*(L1_+2*L2_+2*L3_+L4_));
    end
    %y1
    figure(1)
    plot(t(3000:4000),z1(3000:4000))
    xlabel('时间(s)','fontsize',16)
    ylabel('埀荡量(m)','fontsize',16) 
    title('埀荡量时间曲线','fontsize',16)
    figure(2)
    plot(t(3000:4000),z2(3000:4000))
    xlabel('时间(s)','fontsize',16)
    ylabel('纵摇角(rad)','fontsize',16)
    title('纵摇角时间曲线','fontsize',16)
    

    龙格库塔算法代码引用了两个计算波浪力和波浪扰动力矩的代码Cal_F和Cal_M

    function F_=cal_F(t)
    fi =  pi/3;      % 波浪遭遇角2*pi/3
    B = 1  ;         % width of boat
    d = 0.341 ;      % draft
    delta = 1.76*10^3 ; % displacement
    V = 1.76    ;       % displacement volume
    g = 9.8     ;
    L = 6     ;        % length
    lamda = 6   ;      % 波长
    A_w = 5.448  ;      % area of water plane
    A_m = 0.34   ;     % '中横剖面面积'
    C_w = A_w / (L * B) ; % coefficient of water plane
    C_p = V/(L*A_m) ; %这里是为了保证cp=0.848.所以l乘以0.91,cb也一样的原因
    v = 4.44;
    T = 0.8*sqrt(lamda);
    W_0 = 2*pi/T ;       % 波浪圆频率
    W_e = W_0*(1+v*W_0*cos(fi)/g)  ;  % 波浪遭遇频率
    x_b = -0.23  ;     % 这个很可能是错的,因为我自己变得,按照大尺度船的计算是0.23,但是根据参考文献,这里应该是负的
    M = 1.76*10^3   ;   % weight
    H_0 = B/(2*d);
    a_zz = (0.8*H_0*C_w*M);
    b_zz = ((5.4*(C_w/C_p)*(H_0)^(0.5)-4.7)*delta/((g*L)^0.5));
    a_z0 = (v/W_e^2)*b_zz-a_zz*x_b;
    I_00 = 0.07*C_w*M*L^2;
    a_00 = 0.83*H_0*C_p^2*(0.25*L)^2*(delta/g)^0.5;
    a_0z = -(a_zz*x_b);
    fi = pi/3;
    a = 0.5*0.17*lamda^(3/4);
    k = 2*pi/lamda ;
    n3 = (I_00+a_00)*(M+a_zz)-a_0z*a_z0;
    m1 = 547*cos(5.53*t);
    m2 = 1045*sin(5.53*t);
    F_=(I_00+a_00)/n3*m1-a_z0/n3*m2;
    end
    
    function M_=cal_M(t)
    
    fi =  pi/3;      % 波浪遭遇角2*pi/3
    B = 1  ;         % width of boat
    d = 0.341 ;      % draft
    delta = 1.76*10^3 ; % displacement
    V = 1.76    ;       % displacement volume
    g = 9.8     ;
    L = 6     ;        % length
    lamda = 6   ;      % 波长
    rho = 1.0*10^3;
    A_w = 5.448  ;      % area of water plane
    A_m = 0.34   ;     % '中横剖面面积'
    C_w = A_w / (L * B) ; % coefficient of water plane
    C_p = V/(L*A_m) ; %这里是为了保证cp=0.848.所以l乘以0.91,cb也一样的原因
    C_b = delta/(1000*L*B*d);
    v = 4.44;
    T = 0.8*sqrt(lamda);
    W_0 = 2*pi/T ;       % 波浪圆频率
    W_e = W_0*(1+v*W_0*cos(fi)/g)  ;  % 波浪遭遇频率
    x_b = -0.23  ;     % 这个很可能是错的,因为我自己变得,按照大尺度船的计算是0.23,但是根据参考文献,这里应该是负的
    M = 1.76*10^3   ;   % weight
    H_0 = B/(2*d);
    a_zz = (0.8*H_0*C_w*M);
    b_zz = ((5.4*(C_w/C_p)*(H_0)^(0.5)-4.7)*delta/((g*L)^0.5));
    a_z0 = (v/W_e^2)*b_zz-a_zz*x_b;
    I_00 = 0.07*C_w*M*L^2;
    a_00 = 0.83*H_0*C_p^2*(0.25*L)^2*(delta/g)^0.5;
    a_0z = -(a_zz*x_b);
    fi = pi/3;
    a = 0.5*0.17*lamda^(3/4);
    k = 2*pi/lamda ;
    n1 = (I_00+a_00)-(a_0z*a_z0)/(M+a_zz);
    n2 = rho*g*a*k*exp(-k*d/2)*B*d*sin(k*B*L*sin(fi)/2)/(k*B*L*sin(fi)/2);
    n3 = (I_00+a_00)*(M+a_zz)-a_0z*a_z0;
    m1 = 547*cos(5.53*t);
    m2 = 1045*sin(5.53*t);
    
    M_=(M+a_zz)/n3*m1-a_0z/n3*m2;
    end
    

    Cal_f和Cal_m里面的m1和m2是我用cal_f_m算的,这一步其实就是在算运动方程的波浪干扰力和波浪干扰力矩。

    clc
    syms t x
    fi =  pi/3;      % 波浪遭遇角2*pi/3
    B = 1  ;         % width of boat
    d = 0.341 ;      % draft
    delta = 1.76*10^3 ; % displacement
    V = 1.76    ;       % displacement volume
    g = 9.8     ;
    L = 6     ;        % length
    lamda = 6   ;      % 波长
    rho = 1.0*10^3;
    A_w = 5.448  ;      % area of water plane
    A_m = 0.34   ;     % '中横剖面面积'
    C_w = A_w / (L * B) ; % coefficient of water plane
    C_p = V/(L*A_m) ; %这里是为了保证cp=0.848.所以l乘以0.91,cb也一样的原因
    C_b = delta/(1000*L*B*d);
    v = 4.44;
    T = 0.8*sqrt(lamda);
    W_0 = 2*pi/T ;       % 波浪圆频率
    W_e = W_0*(1+v*W_0*cos(fi)/g)  ;  % 波浪遭遇频率
    x_b = -0.23  ;     % 这个很可能是错的,因为我自己变得,按照大尺度船的计算是0.23,但是根据参考文献,这里应该是负的
    GM_l = (L^2*(5.55*C_w+1)^3/3450)/(C_b*d);
    M = 1.76*10^3   ;   % weight
    H_0 = B/(2*d);
    a_zz = (0.8*H_0*C_w*M);
    b_zz = ((5.4*(C_w/C_p)*(H_0)^(0.5)-4.7)*delta/((g*L)^0.5));
    c_zz = rho*g*A_w ;
    a_z0 = (v/W_e^2)*b_zz-a_zz*x_b;
    b_z0 = -b_zz*x_b+v*a_zz  ;   % 0其实就是theta
    c_z0 = -rho*g*A_w*x_b;
    I_00 = 0.07*C_w*M*L^2;
    a_00 = 0.83*H_0*C_p^2*(0.25*L)^2*(delta/g)^0.5;
    b_00 = 0.08*H_0*delta*L^2/(g*L)^0.5;
    c_00 = rho*g*(V*GM_l);
    a_0z = -(a_zz*x_b);
    b_0z = -(b_zz*x_b);
    c_0z = -rho*g*A_w*x_b-v*b_zz;
    fi = pi/3;
    a = 0.5*0.17*lamda^(3/4);
    k = 2*pi/lamda ;
    n1 = (I_00+a_00)-(a_0z*a_z0)/(M+a_zz);
    n2 = rho*g*a*k*exp(-k*d/2)*B*d*sin(k*B*L*sin(fi)/2)/(k*B*L*sin(fi)/2);
    c_f=n2*int(cos(k*x*cos(fi)-W_e*t),x,-3,3)
    c_m=n2*int(x*cos(k*x*cos(fi)-W_e*t),x,-3,3)
    
    jk1=5728700179944927/(549755813888*pi)
    jk2=3116096361125085/562949953421312
    jk3=17186100539834781/(274877906944*pi^2)
    jk4=3116096361125085/562949953421312
    

    用的时候把这6个代码放在一起,运行龙格库塔算法的那个程序。

    3 看看结果


    在这里插入图片描述

    可以看到都是周期性的运动了,而且符合纵摇-埀荡联合运动的规律。
    一起学习,一起进步!

    展开全文
  • 这篇博客是在将我上一篇博客的matlab代码移植到python中,应为后续要开展深度强化学习下的船舶减摇研究,总的来说还是在python上做这项工作比较合适。具体方程的解法和背景不在赘述,见...

    0 写在前面

    这篇博客是在将我上一篇博客的matlab代码移植到python中,应为后续要开展深度强化学习下的船舶减摇研究,总的来说还是在python上做这项工作比较合适。具体方程的解法和背景不在赘述,见https://blog.csdn.net/Mezikov/article/details/107461970

    1 代码

    import math
    import numpy as np
    import matplotlib.pyplot as plt
    from math import e
    from numpy import reshape
    
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    
    ##################   写在前面:这个代码只适合用于已知浪向角的规则波   ########################
    
    ts, ys, zs, ys_, zs_ = [], [], [], [], []
    # all below is coefficients of wave
    fi = math.pi / 3
    B = 1
    d = 0.341
    delta = 1.76 * 10 ** 3
    V = 1.76
    g = 9.8
    L = 6
    lamda = 1.64
    rho = 1.025 * 10 ** 3
    A_w = 5.448
    A_m = 0.346
    C_w = A_w / (L * B)
    C_p = V / (L * A_m)
    C_b = delta / (1000 * L * B * d)
    v = 4.44
    T = 0.8 * math.sqrt(lamda)
    W_0 = 2 * math.pi / T
    W_e = W_0 * (1 + v * W_0 * np.cos(fi) / g)
    x_b = -0.23
    GM_l = (L ** 2 * (5.55 * C_w + 1) ** 3 / 3450) / (C_b * d)
    M = 1.6 * 10 ** 3
    H_0 = B / (2 * d)
    a_zz = 0.8 * H_0 * C_w * M
    b_zz = (5.4 * (C_w / C_p) * (H_0 ** 0.5) - 4.7) * delta / ((g * L) ** 0.5)
    c_zz = rho * g * A_w
    a_z0 = (v / W_e ** 2) * b_zz - a_zz * x_b
    b_z0 = -b_zz * x_b + v * a_zz
    c_z0 = -rho * g * A_w * x_b
    I_00 = 0.07 * C_w * M * L ** 2
    a_00 = 0.83 * H_0 * C_p ** 2 * (0.25 * L) ** 2 * (delta / g) ** 0.5
    b_00 = 0.08 * H_0 * delta * L ** 2 / (g * L) ** 0.5
    c_00 = rho * g * (V * GM_l)
    a_0z = -(a_zz * x_b)
    b_0z = -(b_zz * x_b)
    c_0z = -rho * g * A_w * x_b - v * b_zz
    fi = math.pi / 3
    a = 0.5 * 0.17 * lamda ** (3 / 4)
    k = 2 * math.pi / lamda
    n1 = (I_00 + a_00) - (a_0z * a_z0) / (M + a_zz)
    n2 = rho * g * a * k * math.e ** (-k * d / 2) * B * d * np.sin(k * B * L * np.sin(fi) / 2) / (
                k * B * L * np.sin(fi) / 2)
    n3 = (I_00 + a_00) * (M + a_zz) - a_0z * a_z0
    
    
    #################   以上全是水动力系数的计算,下面是四阶龙格库塔算法   ######################
    def cal_F(t):
        m1 = 547 * np.cos(5.53 * t)
        m2 = 1045 * np.sin(5.53 * t)
        F_ = (I_00 + a_00) / n3 * m1 - a_z0 / n3 * m2
        return F_
    
    
    def g(y1, z1, y2, z2):
        j = (b_zz * (I_00 + a_00) - a_z0 * b_0z) / n3 * z1 + \
            (c_zz * (I_00 + a_00) - a_z0 * c_0z) / n3 * y1 + \
            (b_z0 * (I_00 + a_00) - a_z0 * b_00) / n3 * z2 + \
            (c_z0 * (I_00 + a_00) - a_z0 * c_00) / n3 * y2
        # 文献结果
        # j= 0.9716*np.cos(1.245*t-0.5566)-0.0674*np.cos(1.245*t-1.4375)-0.5252*z1-1.1424*y1-2.1741*z2-3.5594*y2
        return j
    
    
    def cal_M(t):
        m1 = 547 * np.cos(5.53 * t)
        m2 = 1045 * np.sin(5.53 * t)
        M_ = (M + a_zz) / n3 * m1 - a_0z / n3 * m2
        return M_
    
    
    def g_(y1, z1, y2, z2):
        # m1 = 547*np.cos(5.53*t)
        # m2 = 1045*np.sin(5.53*t)
        # j = (0.2580 * np.sin(5.53 * t) - 0.0042 * np.cos(5.53 * t)) - 0.2137 * z2 - 33.3415 * y2 - 0.0036 * z1 - 0.1641 * y1
        j = (b_00 * (M + a_zz) - a_0z * b_z0) / n3 * z2 + \
            (c_00 * (M + a_zz) - a_0z * c_z0) / n3 * y2 + \
            (b_0z * (M + a_zz) - a_0z * b_zz) / n3 * z1 + \
            (c_0z * (M + a_zz) - a_0z * c_zz) / n3 * y1
        # j=1.1854*np.cos(1.245*t- 1.4375)-0.006*np.cos(1.245*t-0.5566)-7.9154*z2-20.3425*y2-0.0305*z1-0.0547*y1
        return j
    
    def runge_kutta():
        """
            兄弟们,之前这里被网上的资料坑了,我发现按照oringin_try里面的办法,龙格库塔里面只有return的那个量会
            更改,其余变量都不变化,如果这是个一阶微分方程,那只有一个量在变化,可现在是二元二阶方程组啊!我现在改
            过来了,就用这种新方法了。
        """
        t = 0
        dt = 0.1
        y1 = 0
        z1 = 0
        y2 = 0
        z2 = 0
        while t <= 500:
            ys.append(y1)
            zs.append(z1)
            ys_.append(y2)
            zs_.append(z2)
            ts.append(t)
            t += dt
    
            k1 = z1  # f(t,y1,z1,y2,z2)
            k1_ = z2  # f_(t,y1,z1,y2,z2)
            l1 = cal_F(t) - g(y1, z1, y2, z2)
            l1_ = cal_M(t) - g_(y1, z1, y2, z2)
            k2 = z1 + 0.5 * dt * l1
            k2_ = z2 + 0.5 * dt * l1_
            l2 = (cal_F(t) + cal_F(t + dt)) / 2 - g(y1 + 0.5 * dt * k1, z1 + 0.5 * dt * l1, y2 + 0.5 * dt * k1_,
                                                    z2 + 0.5 * dt * l1_)
            l2_ = (cal_M(t) + cal_M(t + dt)) / 2 - g_(y1 + 0.5 * dt * k1, z1 + 0.5 * dt * l1, y2 + 0.5 * dt * k1_,
                                                      z2 + 0.5 * dt * l1_)
            k3 = z1 + 0.5 * dt * l2
            k3_ = z2 + 0.5 * dt * l2_  #  我曹之前这里z2+0.5+dt*l2,我说怎么会这么奇怪一直和Matlab对不上,写代码不能分心!
            l3 = (cal_F(t) + cal_F(t + dt)) / 2 - g(y1 + 0.5 * dt * k2, z1 + 0.5 * dt * l2, y2 + 0.5 * dt * k2_,
                                                    z2 + 0.5 * dt * l2_)
            l3_ = (cal_M(t) + cal_M(t + dt)) / 2 - g_(y1 + 0.5 * dt * k2, z1 + 0.5 * dt * l2, y2 + 0.5 * dt * k2_,
                                                      z2 + 0.5 * dt * l2_)
            k4 = z1 + dt * l3  # wo ta ma zhengde cao le 这里我为什么z1会乘以0.5啊。
            k4_ = z2 + dt * l3_  # f_(t,y1,z1,y2,z2)+0.5*dt*l3_
            l4 = cal_F(t + dt)- g(y1 + dt * k3, z1 + dt * l3, y2 + dt * k3_,
                                                    z2 + dt * l3_)  # 又乘了0.5,第四步应该都是dt而不是0.5*dt
            l4_ = cal_M(t + dt)- g_(y1 + dt * k3, z1 + dt * l3, y2 + dt * k3_, z2 + dt * l3_)
    
            y1 = y1 + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6
            z1 = z1 + dt * (l1 + 2 * l2 + 2 * l3 + l4) / 6
            y2 = y2 + dt * (k1_ + 2 * k2_ + 2 * k3_ + k4_) / 6
            z2 = z2 + dt * (l1_ + 2 * l2_ + 2 * l3_ + l4_) / 6
    
    
    ###################   主程序入口,代码到这里就算编完了  #########################3
    if __name__ == "__main__":
        runge_kutta()
        print(ys)
        # print(ys)
        plt.figure(1)
        plt.plot(ts[2000:3000], zs[2000:3000], label='埀荡')
        plt.xlabel('time(s)')
        plt.ylabel('埀荡量(m)')
        plt.legend()
        plt.show()
        plt.figure(2)
        plt.plot(ts[2000:3000], zs_[2000:3000], label='纵摇')
        plt.xlabel('time(s)')
        plt.ylabel('纵摇角(m)')
        plt.legend()
        plt.show()
        print('最大埀荡量', np.max(ys), )
        print('最大纵摇角', np.max(ys_), )
    

    2 结果

    在这里插入图片描述
    在这里插入图片描述

    3 总结

    但是不知道为什么,matlab算的和python算的结果有一点点(2%左右误差)不一样,这是为什么呢

    展开全文
  • 我们常常用微分方程来描述现实世界中的一些物理现象。由于微分方程的复杂性,只有在很简单的情况下才能得到微分方程的解析。由于计算机的发展,采用数值方法求解微分方程的数值近似得到...微分方程主要分为常微分
  • 对于在反应工程中常见的一类特殊的二阶常微分方程边值问题,给出了二分法初值化求解的一种新方法。具体求解了多孔催化剂和多孔电极两个数学模型,给出了在不同参数下二者的曲线。与传统的打靶法相比,此方法回避了...
  • 对于在反应工程中常见的一类特殊的二阶常微分方程边值问题,给出了二分法初值化求解的一种新方法。具体求解了多孔催化剂和多孔电极两个数学模型,给出了在不同参数下二者的曲线。与传统的打靶法相比,此方法回避了...
  • 大意就是一个力学模型,现在需要解一个二阶常微分方程但是很复杂,大概是有(D2y=f(y)的形式,f是一个由一些向量点乘、叉乘blabla得出来的以y为自变量的函数,但主要是三角函数与多项式、分式函数迭代,或者说,这个...
  • matlab解常微分方程——符号解法

    千次阅读 2019-09-15 10:27:26
    用matlab可以解决许多数学问题,本文章主要讲解利用matlab解常微分方程 1、首先得介绍一下,在matlab中解常微分方程有两种方法,一种是符号解法,另一种是数值解法。在本科阶段的微分数学题,基本上可以通过符号...
  • 练习。使用四阶龙格库塔法求解常微分方程组,通用性较佳。附加一个振动方程求解的案例。振动方程是一个二阶微分方程,转化为两个方程组以后用编写的代码求解。
  • 本文主要介绍matlab中求解常微分方程(组)的dsolve和ode系列函数,并通过例子加深读者的理解。 一、符号介绍  D: 微分符号;D2表示二阶微分,D3表示三阶微分,以此类推。 二、函数功能介绍及例程 1、dsolve ...
  • MATLAB解决常微分方程

    2018-08-14 11:38:00
    首先得介绍一下,在matlab解常微分方程有两种方法,一种是符号解法,另一种是数值解法。在本科阶段的微分数学题,基本上可以通过符号解法解决。 用matlab解决常微分问题的符号解法的关键命令是dslove...
  • Matlab学习——求解微分方程(组)

    万次阅读 多人点赞 2018-05-29 17:18:58
    函数 dsolve 用来解决常微分方程(组)的求解问题,调用格式为 X=dsolve(‘eqn1’,’eqn2’,…) 如果没有初始条件,则求出通,如果有初始条件,则求出特 系统缺省的自变量为 t。 2.函数d...
  • 常微分数值解matlab代码symODE2:具有多项式系数的二阶常微分方程的符号分析 作者:Tolga Birkandan 电子邮件:伊斯坦布尔技术大学物理系,土耳其伊斯坦布尔 34469。 详情请参阅。 提出了一种用于对具有多项式系数的...
  • 本文主要介绍matlab中求解常微分方程(组)的dsolve和ode系列函数,并通过例子加深读者的理解。 一、符号介绍    D: 微分符号;D2表示二阶微分,D3表示三阶微分,以此类推。 二、函数功能介绍及例程 1、...
  • matlab解微分方程

    万次阅读 多人点赞 2018-02-07 16:07:01
    一般来说,在matlab解常微分方程有两种方法,一种是符号解法,另一种是数值解法。在本科阶段的微分数学题,基本上可以通过符号解法解决。 用matlab解决常微分问题的符号解法的关键命令是dslove命令。该命令中可以...
  • 首先用高数知识求解非齐次系数微分方程 再利用信号与系统中冲激响应求解验证 利用MATLAB求解验证 y=dsolve('D2y+3*Dy+2*y=exp(-t)','y(0)=1','Dy(0)=2','t') 得出结果: y = (t - 2 exp(-t) + 3) exp(-t) ...
  • 表1 解常微分方程主要MATLAB指令 主题词 意义 主题词 意义 ode45 4、5阶Runge-kutta法 ode23s 刚性方程组二阶Rosenbrock法 ode23 2、3阶Runge-kutta法 ode23tb 刚性方程组低精度算法 ...
  • 运用追赶法来三对角线性方程组...用差分法求解二阶常微分方程边值问题时,最后常规为求解具有三对角系数矩阵的线性方程组.对三对角矩阵实行Doolittle(或Crout)分解,便得到求解三对角方程组的最有效方法---追赶法.

空空如也

空空如也

1 2 3
收藏数 46
精华内容 18
关键字:

matlab解二阶常微分方程

matlab 订阅