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

3 次查看（过去 30 天）
Satish Jawalageri
Satish Jawalageri 2020-7-20
评论： Satish Jawalageri 2020-7-22
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')
  3 个评论
显示 1更早的评论
Satish Jawalageri
Satish Jawalageri 2020-7-20
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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Serhii Tetora
Serhii Tetora 2020-7-22
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by