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 ont = 0:0.1:500for 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)];endprofile offprofile 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')
Best Answer