Tuesday, March 21, 2017

                                                              Chapter-10
    SOLUTION OF DIFFERENTIAL EQUATION USING MATLAB
MATLAB has large library of tools that can be used for solving differential equation.
dy/dt = f(t,y)  for t0<t<t with y = y0   at  t = t0
We can solve ordinary differential equation using MATLAB programming. We should note that using MATLAB we cannot find y(t) (or y as a function of t). We can find values of y at specified values of t.
There are following steps for solving differential equation:
  • ·         First step is to create a user-defined function file that calculates dy/dt for given values of t and y.

Syntax of function file:
function dydt = name_of_function_file( input_args )
dydt = f(t,y);
end
  • ·         Second step is to select method of solution.

Select the numerical method that you would like MATLAB to use in the solution. Several of the numerical methods are available as built-in functions in MATLAB. There are following MATLAB built-in functions that are used for solving differential equation:
1.       ode45 : It is one step solver and based on explicit Runge-kutta method.
2.       ode23 : It is also one step solver and based on explicit Runge-kutta method. It is often quicker but less accurate than ode45.
  • ·         Third step is to solve the differential equation using MATLAB.

Syntax for using command that is used to solve an initial value ordinary differential equation is:
[t,y] = built-in_function('name_of_function_file',tspan,y0)



For example:
We have differential equation:
dy/dt = (t3 – 2y)/t  for 1<t<3 with y = 4.2 at t = 1
Function file for following differential equation:
function dydt = ODE1(t,y)
%function for solving differential eqquation
%   Detailed explanation goes here
dydt = (t^3 - 2*y)/t;
end

Suppose We have selected ode45 method for solving differential equation. Script file for above differential equation:
%function for solving differential equation
[t y] = ode45('ODE1’,[1:0.01:3],4.2);%we have got values of y and t in workspace after the execution of this command.
plot(t,y)%for plotting graph of y with respect to t

Q) Solve the following differential equation using MATLAB:

dy/dx  =  (x2 - 2x + 3)/(y2)  for  0.5<x<3  with y(0.5)  =  2

function dydx = ODE3( x,y )
%fuunction file for solving differential equation
dydx = (x^2 - 2*x + 3)/(y^2);
end

clc
clear all
close all
xspan = [0.5:0.01:3];
y1 = 2;
[x y] = ode45('ODE3',xspan,y1);
plot(x,y)




                                       APPLICATION IN CALCULUS

An equation with one variable can be written in the form f(x) = 0. The solution is the value of x where the function crosses the x-axis (the value of function is zero), which means that the function changes sign at x. Using MATLAB, we can easily find the solution of function.
Syntax for finding solution:

‘function’ is the function to be solved. It can be entered in three different ways:
  • ·         The simplest way is to enter the mathematical expression as string.
  • ·         The function can also be created as a user-defined function in a function file, then name of the function typed as a string.

X0 can be a scalar or a two-element vector. If it is entered as a scalar. It has to be a point near which the function crosses the x-axis. If x0 is entered as a vector of two elements, function should have different sign at these two elements.
Q) Using MATLAB, determine the solution of equation?
                                  x*exp(-x) – 0.2 = 0

fplot('x*exp(-x) - 0.2',[0 8])

By executing above command in MATLAB, we can see the graph of function and from graph we conclude that function has one solution near 0.7 and another solution near 2.8.
>> x1 = fzero('x*exp(-x) - 0.2',0.7)
x1 =

    0.2592
>> x2 = fzero('x*exp(-x) - 0.2',2.8)
x2 =

    2.5426
Therefore, by executing above two commands in MATLAB, we have concluded that function has one solution at 0.2592 and another at 2.5426.
  • ·         The fzero command finds zeroes of a function only where the function crosses the x-axis. It does not find zero at points where the function touches but does not cross the x-axis.
  • ·         [x fval] = fzero('function',x0) assigns the value of the function at x to the variable fval.

Q) Determine the first three positive roots of equation:    2*sin(x) - x.^0.5 + 2.5 = 0 by using MATLAB?

clc
close all
clear all
%By executing following graphs in matlab and by watching following graphs
%we will get approximate ideas about roots
%fplot('2*sin(x)-x.^0.5+2.5',[3 4]);
%fplot('2*sin(x)-x.^0.5+2.5',[6 7]);
%fplot('2*sin(x)-x.^0.5+2.5',[9 10]);
%from first above graph we conclude that one root is near 3.5
%from second above graph we conclude that one root is near 6.3
%from third graph we conclude that one root is near 9.15
x1 = fzero('2*sin(x)-x.^0.5+2.5',3.5);
x2 = fzero('2*sin(x)-x.^0.5+2.5',6.3);
x3 = fzero('2*sin(x)-x.^0.5+2.5',9.15);
fprintf('function has solution at points %f,%f and %f',x1,x2,x3);


*link for downloading program threepositiveroot.m from folder calculus
FINDING A MINIMUM OR A MAXIMUM OF A FUNCTION
In MATLAB, the value of x where a one-variable function f(x) within the interval x1<x<x2 has a minimum can be determined with the fminbnd command which has the form:




     









      
  •       The function can be entered as a string, as the name of a function file, or as the name of an inline function, in the same way as with the fzero command.
  • ·         The value of the function at the minimum can be added to the output by using the option:


[x fval] = fminbnd('function',x1,x2)
Where the value of the function at X is assigned to the variable fval
Consider the function f(x) = x^3 -12*x^2 + 40.25*x -36.5

fplot('x^3 - 12*x^2 + 40.25*x - 36.5',[0 8])

By executing above command in command window, we can see the graph of above function.

>> [x fval] = fminbnd('x^3 - 12*x^2 + 40.25*x - 36.5',0,8)

x =
  5.6073
fval =

  -11.8043

Therefore minimum of above function at x = 5.6073 and value of function at this point is -11.8043.
We should note that fminbnd command gives local minima and not the absolute minimum which is at x=0.
  • ·         fminbnd command tries to find local minima and not the absolute minima. If a local minimum is not found, the function returns the end point (x1 or x2) with the lower function value.
  • ·         The fminbnd command can also be used to find the maximum of a function. This is done by multiplying the function by -1 and finding the minimum.

NUMERICAL INTEGRATION
 Integration is a common mathematical operation in Science and Engineering. It is assumed that reader has basic knowledge of integration.
There are two MATLAB built-in functions, quad and trapz that are used for finding integration. Syntax for using quad command:


The function can be entered as string, as the name of a function file in the same way as with the fzero command.
Q) Use MATLAB to calculate the integration of function f(x) = 1/(0.8x2 + 0.5x+2) from x = 0 to x = 8?
                                                                            OR
Q) Use MATLAB to calculate the area bounded by curve y = 1/(0.8x2 + 0.5x+2), x=0, x=8 and x-axis?

clc
clear all
a = quad('(0.8*x.^2+0.5*x+2).^-1',0,1);

*link for downloading program iintegration.m from folder intdiffeq
Therefore by using quad command we can find integration in just one step.

Trapz command:
The trapz command can be used for integrating a function that is given as data points. It uses numerical trapezoidal method of integration. Syntax for using command:
trapz(x,y)
where x and y are vectors with the x and y coordinates of the points respectively. The vectors must be of same length. For example:

>> X = [0 1 2 3 4 5];
>> Y = [0 1 2 3 4 5];
>> trapz(X,Y)

ans =

   12.5

Therefore area bounded by curve X-Y is equal to 12.5. We can easily verify this from the X-Y plot shown below.


























                                                    Chapter-8
                                       THREE-DIMENSIONAL PLOTS
A three-dimensional line plot is a line that is obtained by connecting points in a three-dimensional space. A basic 3-D plot is created by plot3 command. Syntax for using plot3 command:


Line color specifier is an optional property that can be used to define the color of line. The line color specifiers are following:

Line color
Specifier
Red
Green
Blue
Cyan
Magenta
Blue
Black
yellow
r
g
b
c
m
b
k
y

LineWidth is an optional property that specifies the width of line. Default value of LineWidth is 0.5 and we can specify its value in the command. The three vectors x, y and z must have the same number of elements.
If x, y and z are given as function of parameter t:

x = (2 + 4*cos(t)).*cos(t);
y = (2 + 4*cos(t)).*sin(t);
z = t.^2;

The xlabel and ylabel commands:
Labels can be placed next to the axis with the xlabel and ylabel commands which have the form:
xlabel('text as string');
ylabel('text as string');
zlabel('text as string');
The grid command:
grid on     Add grid lines to the plot
grid off     Removes grid lines from the plot
A plot of points for t = 0:0.01:20 can be produced by following script file:

%3D graph
%use of plot3 command to plot the graph
clc
clear all
t = 0:0.01:20;
x = (2 + 4*cos(t)).*cos(t);
y = (2 + 4*cos(t)).*sin(t);
z = t.^2;
plot3(x,y,z,'k','linewidth',1)
grid on%add grid lines to the plot
xlabel('x');
ylabel('y');
zlabel('z');

*link for downloading program threedgraph.m from folder 3Dgraph


By running following script file in MATLAB we will get idea on 3-D plotting in MATLAB using plot3 command.
MESH AND SURFACE PLOTS:
Mesh and surface plots are three-dimensional plots that are used for plotting functions of form z = f (x, y) where the x and y are the independent variables and z is the dependent variable.
Mesh and surface plots are created in three steps:
  • ·         Divide the domain of independent variables x and y.
  • ·         form the matrices of independent variables.
  • ·         Calculate the value of dependent variable z at each point in domain.
  • ·         Create the plot.

MATLAB has in-built command called meshgrid, that can be used for creating X and Y matrices. Syntax for using meshgrid command:



Fig: 8.1 shows how to create the matrices of independent variables
>> x = -1:4;
>> y = 1:5;
>> [X,Y] = meshgrid(x,y)
X =

    -1     0     1     2     3     4
    -1     0     1     2     3     4
    -1     0     1     2     3     4
    -1     0     1     2     3     4
    -1     0     1     2     3     4
Y =

     1     1     1     1     1     1
     2     2     2     2     2     2
     3     3     3     3     3     3
     4     4     4     4     4     4
     5     5     5     5     5     5

The value of z at each point can be created by using element-by-element calculations in the same way it is used with vectors.
Mesh or surface plot is created with mesh or surf commands, which have the form:
surf(X,Y,Z)                    mesh(X,Y,Z)
Where X and Y are matrices of independent variable and Z is a matrix with value at each point in the domain of x and y.
As an example, the following script file contains a complete program that plots the function:


%use of mesh command and use of surf command
x = -3:0.25:3;
y = -3:0.25:3;
[X,Y] = meshgrid(x,y);
Z = 1.8.^(-1.5*sqrt(X.^2+Y.^2)).*sin(X).*cos(0.5*Y);
subplot(1,2,1),surf(X,Y,Z)
subplot(1,2,2),mesh(X,Y,Z)
%by watching the output we can clearly saw the difference between surf and
%mesh command and subplot command is used for plotting two graph on same %sheet
xlabel('x')
ylabel('y')
zlabel('z')


*link for downloading program: meshcommand.m from folder 3Dgraph
OTHER COMMANDS USED FOR THREE DIMENSIONAL PLOTTING:
meshz and meshc command: meshz command draws curtain around the mesh and meshc commands draws a contour plot beneath the mesh.



Q) plot the function
over the domain -3<x<3 and -3<y<3 using meshz and meshc command and observe the difference between two commands?



%use of meshz command and use of meshc command
x = -3:0.25:3;
y = -3:0.25:3;
[X,Y] = meshgrid(x,y);
Z = 1.8.^(-1.5*sqrt(X.^2+Y.^2)).*sin(X).*cos(0.5*Y);
subplot(1,2,1),meshc(X,Y,Z)%draws a contour plot beneath the mesh
subplot(1,2,2),meshz(X,Y,Z)%draws a curtain around the mesh
%by watching the output we can clearly saw the difference between meshz and
%meshc command
xlabel('x')
ylabel('y')
zlabel('z')


*link for downloading program meshzandmeshc.m


surfc and surfl command: surfc command draws a contour plot beneath the surface and surfl command draws surface plot with lighting






Q)  plot the function

over the domain -3<x<3 and -3<y<3 using surfc and surfl command and observe the difference between two commands?



%comparison between surf and surfc command
x = -3:0.25:3;
y = -3:0.25:3;
[X,Y] = meshgrid(x,y);
Z = 1.8.^(-1.5*sqrt(X.^2+Y.^2)).*sin(X).*cos(0.5*Y);
subplot(1,2,1),surfc(X,Y,Z)%draws a contour plot beneath the surface
subplot(1,2,2),surfl(X,Y,Z)
%by watching the output we can clearly saw the difference between surfc and
%surf command
xlabel('x')
ylabel('y')
zlabel('z')
 *link for downloading program surflandsurfc.m from folder 3Dgraph


The view command:

The view command controls the direction from which the plot is viewed. This is done by specifying the direction in terms of azimuthal and elevation angles or by defining points in space from which the plot is viewed. view command has syntax:
view(azimuthal angle,elevation angle) or view([azimuthal angle,elevation angle])
  • ·         Azimuthal angle in the x-y plane is measured relative to the negative y axis direction and defined as positive in counter-clockwise direction.
  • ·         Elevation angle is measured from the x-y plane. A positive value corresponds to opening an angle in the direction of z axis.
  • ·         The default azimuthal angle = -37.5 and default elevation angle = 30.

By choosing appropriate azimuthal and elevation angle the view command can be used to plot projections of 3-D plots on various planes according to following table:
Projection plane
az value
el value
x-y(top view)
0
90
x-z(side view)
0
0
y-z(side view)
90
0

Q) Write three different programs to Draw the projection of function 

 on x-y plane, y-z plane and x-z plane?
First program: projection on x-y plane


clc
clear all
close all
x = -3:0.25:3;
y = -3:0.25:3;
[X,Y] = meshgrid(x,y);
Z = 1.8.^(-1.5*sqrt(X.^2+Y.^2)).*sin(X).*cos(0.5*Y);
surf(X,Y,Z)
view(0,90)
 *link for downloading program projectiononxy.m from folder 3D graph




Second program: projection on y-z plane
%projection on yz plane of function by using view command
clc
clear all
close all
x = -3:0.25:3;
y = -3:0.25:3;
[X,Y] = meshgrid(x,y);
Z = 1.8.^(-1.5*sqrt(X.^2+Y.^2)).*sin(X).*cos(0.5*Y);
surf(X,Y,Z)
view(90,0)

*link for downloading program projectiononyz.m from folder 3D graph




Third program: projection on x-z plane
%projection on xz plane of function by using view command
clc
clear all
close all
x = -3:0.25:3;
y = -3:0.25:3;
[X,Y] = meshgrid(x,y);
Z = 1.8.^(-1.5*sqrt(X.^2+Y.^2)).*sin(X).*cos(0.5*Y);
surf(X,Y,Z)
view(0,0)

*link for downloading program projectiononxzplane.m from folder 3D graph



Q) four charges of magnitude 2×10-10 C are placed at point (0.25,0,0) (-0.25,0,0) (0,0.25,0) and                    (0,-0.25,0). plot the electrical potential due to the particles at points in x-y plane that are located in domain -0.5<x<0.5 -0.5<y<0.5. Make the plot such that the x-y plane is the plane of the points and the
z-axis is the magnitude of electrical potential?


clc
clear all
eps0 = 8.85e-12;
q1 = 2e-10;
q2 = 2e-10;
q3 = 2e-10;
q4 = 2e-10;
k = 1/(4*pi*eps0);
x = -0.5:0.01:0.5;
y = -0.5:0.01:0.5;
[X,Y] = meshgrid(x,y);
r1 = sqrt((X+0.25).^2 + Y.^2);
r2 = sqrt((X-0.25).^2 + Y.^2);
r3 = sqrt((Y+0.25).^2 + X.^2);
r4 = sqrt((Y-0.25).^2 + X.^2);
V = k*(q1./r1 + q2./r2 + q3./r3 + q4./r4);
mesh(X,Y,V)
xlabel('x');
ylabel('y');
zlabel('z');
  

*link for downloading program potential.m from folder 3Dgraph