ENGG952 Engineering Computing-Implicit Method 

Mass-Spring-Damper system for an Ordinary Differential Equation Free vibrations in a mechanical system are caused by initial conditions of the parameters: displacement, velocity, or acceleration where no external force is allowed to interact. The mechanical system in the free vibration will oscillate with its natural frequency and eventually settle down to zero as a result of damping effects. 

Dividing through by mass, m,The natural frequency of the system and the damping ratio are replaced as shown in the following equation,The solution to the differential equation. 

The transient solution of the system is given as,The constants A and B are obtained when the initial conditions are inserted. The final solution is obtained as,From the initial conditions, 

The Runge-Kutta Heun Method considering a certain range or period of time

tspan=[0:0.1:50];

y0=[0.02;0]; 

[t,y]=ode45('unforced2',tspan,y0);

plot(t,y(:,1));

grid on

xlabel('time')

ylabel('Displacement')

title('Critically-damped Modelling of MSD system (C=40)') 

%% order 647236

tspan=[0:0.1:50];

y0=[0.02;0];

[t,y]=ode45('unforced1',tspan,y0);

plot(t,y(:,1));

grid on

xlabel('time')

ylabel('Displacement')

title('Under-damped Modelling of MSD system (C=5)')

hold on;

%% order 647236

function yp=unforced2(t,y)

c=40;

m=20;

k=20;

yp= [y(2); (-((c/m)*y(2))-((k/m)*y(1)))];

%% order 647236

function yp=unforced1(t,y)

c=5;

m=20;

k=20;

yp= [y(2); (-((c/m)*y(2))-((k/m)*y(1)))];

The graphical illustration of the behavior of the mass-spring-damper for an under-damped and critically-damped system are shown in the same axes. The critically damped reduces to zero faster than the under-damped. The damping coefficient determines the rate at which the signal reduces its amplitude or magnitude as it travels. heat conduction in an insulated rod. The non-dimensional form is expressed.

The explicit finite difference method and the implicit crank-Nicholson method are used to solve this equation. The explicit method of parabolic PDE is simple to program though governed by a criterion for the solution to converge. The steps used are restricted by the criterion. There are implicit methods employed as well as they are more stable. There are two implicit method namely the simple implicit method and the Crank-Nicholson method which are simple with low accuracy and complex with high accuracy respectively. The explicit schemes are conditionally stable for the parabolic partial differential equations. The case study is on heat conduction equation over time variable with one spatial dimension.

• Implicit method 

The elliptic parabolic PDEs can be solved using the finite Difference Method. The dimension of time in this case is open and the initial condition is set at t=0.

Step 1: the temperature for the next time step is calculated and indexed as shown below,

For the i^th node: -

Space domain at time l: Time domain at node i:

The time derivative, unfortunately, is less accurate than the spatial derivative. The central difference in terms of x, indexed by I,To ensure stability the time step is limited by a space interval such that, 

• Crank-Nicholson method 

Clear all; % clear memory of variables and figures

alfa = 0.835;

x0 = 0;                             % starting x

xm = 10;                            % ending x

t0 = 0;                             % starting t

tn = 0.5;                            % ending t

m = 5;                              % number of interval in x direction

n = 300;                            % number of interval in t direction

x = linspace(x0,xm,m+1);            % x coordinates of the nodes

t = linspace(t0,tn,n+1);            % t coordinates of the nodes

del_x = (xm - x0)/m;                % delta x

del_t = (tn - t0)/n;                % delta t

hh = del_x*del_x;                   % constant

r = hh/del_t;                       % constant

p = 2*r/alfa;                       % constant

%%%%%   Setting the initial conditions    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%f = inline('sin(x)');              % initial condition expression

T(1,1:m+1) = 0;                     % set T at t = 0

%%%%%   Setting the boundary conditions   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

T(2:n+1,1) = 100.0;                 % T at x = x0

T(2:n+1,m+1) = 50.0;                % T at x = xm

%%%%% Setting the matrix form and solve them at each time step

A_left(1:m-1,1:m-1) = 0;

A_right(1:m-1,1:m-1) = 0;

for i = 1:m-1

    A_left(i,i) = -2 - p;

    A_right(i,i) = 2 - p;

    if i < m-1;

        A_left(i,i+1) = 1;

        A_right(i,i+1) = -1;

    end

    if i > 1;

        A_left(i,i-1) = 1;

        A_right(i,i-1) = -1; 

for j = 1:n

    b = A_right*T(j,2:m)';

    b(1) = b(1) - T(j,1) - T(j+1,1);

    b(m-1) = b(m-1) - T(j,m+1) - T(j+1,m+1);

    T(j+1,2:m) = A_leftb;

%this figure is only suitalbe for dt=0.1sec only.

figure(1)

x=x0:del_x:xm

plot(x,T(11,:),x,T(31,:), x,T(51,:), x,T(101,:),x,T(201,:),x,T(301,:))

legend('t=0.0s','t=0.1s','t=0.2s','t=0.3s','t=0.4s','t=0.5s');

xlabel('x-direction (Non-dimensional x)'); ylabel('temperature (non dimensional temperature)') 

To calculate the temperatures for the first step seeks to move from the known to the unknown,