Taking too long to get the output (Wilson theta method)

Could I know where I might be wrong in the following code as it is taking too long get the result.

WilsonMethod.m

function [depl,vel,accl,t] = WilsonMethod(M,K,C,R)

clc ;
sdof = length(K) ;
% Time step and time duration
ti = 0. ;
dt = 0.1 ;  
tf = 500 ;
t = ti:dt:tf ;
nt = length(t) ;

% Initialize the displacement,velocity and acceleration matrices
depl = zeros(sdof,nt) ;
vel = zeros(sdof,nt) ;
accl = zeros(sdof,nt) ;
Reff = zeros(sdof,nt) ;

% Initial conditions
depl(:,1) = zeros ;
vel(:,1) = zeros ;
accl(:,1) =zeros;% M\(R-K*depl(:,1)-C*vel(:,1)) ;

% Integration constants
tita = 1.4 ;     % Can be changed
a0 = 6/(tita*dt)^2 ; a1 = 3/(tita*dt) ; a2 = 2*a1 ;
a3 = tita*dt/2 ;     a4 = a0/tita ;       a5 = -a2/tita ;
a6 = 1-3/tita ;      a7 = dt/2 ;   a8 = dt^2/6 ;

% Form Effective Stiffness Matrix
Keff = K+a0*M+a1*C ;

%Time step starts
for it = 1:nt-1
    % Calculating Effective Load
    Reff(:,it) = R(:,it)+tita*(R(:,it)-R(:,it))+M*(a0*depl(:,it)+a2*vel(:,it)+2*accl(:,it))+....
                C*(a1*depl(:,it)+2*vel(:,it)+a3*accl(:,it)) ;
    % Solving for displacements at time (t+dt)
    depl(:,it+1) = Keff\Reff(:,it) ;
    % Calculating displacements, velocities and accelerations at time t+dt
    accl(:,it+1) = a4*(depl(:,it+1)-depl(:,it))+a5*(vel(:,it))+a6*accl(:,it) ;
    vel(:,it+1) = vel(:,it)+a7*(accl(:,it+1)+accl(:,it)) ;
    depl(:,it+1) = depl(:,it)+dt*vel(:,it)+a8*(accl(:,it+1)+2*accl(:,it)) ;
end

Ex.m

MA = [8070000,0,-629468070;0,8070000,112980;-629468070,112980,6.800000000000000e+10];
Ad = [8.053095218400001e+06,0,-4.831857131040000e+08;0,2.167940435676214e+05,0;-4.831857131040000e+08,0,3.865485704832000e+10];
Ca = [0,0,0;0,3.241885080000000e+05,0;0,0,1.301151158327999e+09];
Cm = [4.12e+04,0,-2.82e+06;0,1.19e+04,0;-2.82e6,0,3.11e+08];
M = MA+Ad;
K = Ca+Cm;
C = zeros(size(K)) ;        % Damping Matrix

Fg = -79086000;                 %Gravitational force
Fbuoy = 7.844814740000000e+07;  %Buoyancy force
Fp = 2.712318560000001e+06;     %Heave force
profile on
t = 0:0.1:500
for i = 1:length(t)
    Fh = hydro(t(i))
    FhT = transpose(Fh)
    R(:,i) = [-334731.8545+27939.6+6.5*10^5+FhT(i);-3517000+Fg+Fbuoy+Fp;-112510430.2+3.44*10^6+266.5*10^5+(FhT(i)*18)];
end
profile off
profile viewer
[depl,vel,accl,t] = WilsonMethod(M,K,C,R) ;
depl' 
figure(1), clf 
plot(t,depl(1,:)), xlabel('time(s)'), ylabel('surge(m)')
title ('Surge vs Time')

figure(2), clf 
plot(t,depl(2,:)), xlabel('time(s)'), ylabel('heave(m)')
title ('heave vs Time')

figure(3), clf 
plot(t,depl(3,:)), xlabel('time(s)'), ylabel('Pitch(deg)')
title ('Pitch vs Time')

[Serhii Tetora]

Please attach function hydro

[Satish Jawalageri]

Thanks for your reply. Here is the code:

hydro.m

function [Fh] = hydro(t)

Cd = 0.6;
Ca = 0.9699;
R = 4.7;

t=0:0.1:500;

for i=1:length(t)
    vel = compute_wavevelocity(t(i));
    acc = compute_waveacceleration(t(i));
    fun = sym((0.5*997*Cd*2*R*abs(vel)*vel)+(Ca*997*pi*(R^2)*acc)+(pi*(R^2)*997*acc));
    Fh(1,i) = double(int(fun,0,-120));
    
end
plot(t,Fh), xlabel('time(s)'), ylabel('Fh')
title ('Fh vs Time')

end

compute_wavevelocity.m:

function velocity = compute_wavevelocity(t)
       T = 10;
       WN = (9.8*(T^2))/(2*pi());
       k = (2*pi())/WN;
       omega = (2*pi())/T;
       y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
       d = 320;
       H = 4;
       vel = (omega*H/2)* ((cosh(k*y))/(sinh(k*d)))* cos((k*0)-(omega*t));
       
       velocity = mean(vel,1);
       
end

compute_waveacceleration.m:

function acceleration = compute_waveacceleration(t)
        T = 10;
       WN = (9.8*(T^2))/(2*pi());
       k = (2*pi())/WN;
       omega = (2*pi())/T;
       y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
       d = 320;
       H = 4;
        acc = ((omega^2)*H/2)* ((cosh(k*y))/(sinh(k*d)))* sin((k*0)-(omega*t));
        
        acceleration = mean(acc,1);
               
    end

[Satish Jawalageri]

Any suggestions?

Answer 1

[Serhii Tetora]

Without loops it runs faster. Please check it carefully

clc;close all;clear;
MA = [8070000,0,-629468070;0,8070000,112980;-629468070,112980,6.800000000000000e+10];
Ad = [8.053095218400001e+06,0,-4.831857131040000e+08;0,2.167940435676214e+05,0;-4.831857131040000e+08,0,3.865485704832000e+10];
Ca = [0,0,0;0,3.241885080000000e+05,0;0,0,1.301151158327999e+09];
Cm = [4.12e+04,0,-2.82e+06;0,1.19e+04,0;-2.82e6,0,3.11e+08];
M = MA+Ad;
K = Ca+Cm;
C = zeros(size(K));            % Damping Matrix

Fg = -79086000;                 %Gravitational force
Fbuoy = 7.844814740000000e+07;  %Buoyancy force
Fp = 2.712318560000001e+06;     %Heave force

t = 0:0.1:500;
Fh = hydro(t);
FhT = transpose(Fh);

R(1,:) = -334731.8545 + 27939.6 + 6.5e5 + FhT;
R(2,:) = -3517000 + Fg + Fbuoy + Fp;
R(3,:) = -112510430.2 + 3.44e6 + 266.5e5 + 18*FhT;

[depl,vel,accl,t] = WilsonMethod(M,K,C,R) ;
% depl';
figure(1), clf
plot(t,depl(1,:)), xlabel('time(s)'), ylabel('surge(m)'), grid on
title ('Surge vs Time')

figure(2), clf
plot(t,depl(2,:)), xlabel('time(s)'), ylabel('heave(m)'), grid on
title ('heave vs Time')

figure(3), clf
plot(t,depl(3,:)), xlabel('time(s)'), ylabel('Pitch(deg)'), grid on
title ('Pitch vs Time')

function [depl,vel,accl,t] = WilsonMethod(M,K,C,R)

sdof = length(K) ;
% Time step and time duration
ti = 0. ;
dt = 0.1 ;
tf = 500 ;
t = ti:dt:tf ;
nt = length(t) ;

% Initialize the displacement,velocity and acceleration matrices
depl = zeros(sdof,nt) ;
vel = zeros(sdof,nt) ;
accl = zeros(sdof,nt) ;
Reff = zeros(sdof,nt) ;

% Initial conditions
depl(:,1) = zeros ;
vel(:,1) = zeros ;
accl(:,1) =zeros;% M\(R-K*depl(:,1)-C*vel(:,1)) ;

% Integration constants
tita = 1.4 ;     % Can be changed
a0 = 6/(tita*dt)^2 ; a1 = 3/(tita*dt) ; a2 = 2*a1 ;
a3 = tita*dt/2 ;     a4 = a0/tita ;       a5 = -a2/tita ;
a6 = 1-3/tita ;      a7 = dt/2 ;   a8 = dt^2/6 ;

% Form Effective Stiffness Matrix
Keff = K+a0*M+a1*C ;

%Time step starts
for it = 1:nt-1
    % Calculating Effective Load
    Reff(:,it) = R(:,it)+tita*(R(:,it)-R(:,it))+M*(a0*depl(:,it)+a2*vel(:,it)+2*accl(:,it))+....
        C*(a1*depl(:,it)+2*vel(:,it)+a3*accl(:,it)) ;
    % Solving for displacements at time (t+dt)
    depl(:,it+1) = Keff\Reff(:,it) ;
    % Calculating displacements, velocities and accelerations at time t+dt
    accl(:,it+1) = a4*(depl(:,it+1)-depl(:,it))+a5*(vel(:,it))+a6*accl(:,it) ;
    vel(:,it+1) = vel(:,it)+a7*(accl(:,it+1)+accl(:,it)) ;
    depl(:,it+1) = depl(:,it)+dt*vel(:,it)+a8*(accl(:,it+1)+2*accl(:,it)) ;
end
end

function [Fh] = hydro(t)
Cd = 0.6;
Ca = 0.9699;
R = 4.7;

vel = compute_wavevelocity(t);
acc = compute_waveacceleration(t);
fun = sym((0.5*997*Cd*2*R*abs(vel).*vel)+(Ca*997*pi*(R^2).*acc)+(pi*(R^2)*997.*acc));
Fh = double(int(fun,0,-120));
end

function velocity = compute_wavevelocity(t)
T = 10;
WN = (9.8*(T^2))/(2*pi);
k = 2*pi/WN;
omega = 2*pi/T;
y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
d = 320;
H = 4;
vel = (omega*H/2)* ((cosh(k*y))/(sinh(k*d)))* cos((k*0)-(omega.*t));

velocity = mean(vel,1);
end

function acceleration = compute_waveacceleration(t)
T = 10;
WN = (9.8*(T^2))/(2*pi);
k = 2*pi/WN;
omega = 2*pi/T;
y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
d = 320;
H = 4;
acc = ((omega^2)*H/2)* ((cosh(k*y))/(sinh(k*d)))* sin((k*0)-(omega.*t));

acceleration = mean(acc,1);
end

[Satish Jawalageri]

I got it. Thanks a lot.