Kishor Vaigyanik Protsahan Yojana (Department of Science and Technology, Government of India) June 11, 2009

Summer Camp Report on

Numerical Solution to Ordinary Dierential Equations in Scilab

by

Rahul Kumar Soni ([email protected])

Department of Fuel and Mineral Engineering Indian School of Mines University Dhanbad under the mentorship of

Professor Kannan M. Moudgalya ([email protected])

Department of Chemical Engineering Indian Institute of Technology Bombay, India 1

Certicate This is to certify that Mr. Rahul Kumar Soni, a student of Department of Fuel and Mineral Engineering from Indian School of Mines University Dhanbad has done his Kishor Vaigyanik Protsahan Yojana summer camp project at Department of Chemical Engineering of Indian Institute of Technology Bombay from May 11, 2009 to June 10, 2009 under my guidance. The project work entitled Numerical Solution to Ordinary Dierential Equations in Scilab embodies the original work done by Mr. Rahul Kumar Soni during his summer camp.

Signature Professor Kannan M. Moudgalya Department of Chemical Engineering Indian Institute of Technology Bombay

2

Acknowledgement My sincere thanks to Prof. Kannan Moudgalya for having given me this opportunity to work on such a lucrative project. I wish to thank the KVPY cell for having provided a great work environment. This project would not have been complete without the guidance and timely inputs from Prof. G K Srinivasan ( Professor, Department of Mathematics, IIT Bombay), Prof. V K Gupta (Tata Steel Chair Professor at Department of Fuel and Mineral Engineering, ISMU), Dr. Nikkam Suresh(Asso. Prof. and HOD, Department of Fuel and Mineral Engineering, ISMU), Dr. Biswajit Paul (Department of Environment Science and Engineering, ISMU), Dr.

S Bhattacharya

(Asso. Prof., Department of Fuel and Mineral Engineering, ISMU), Dr.

N K Singh(Asso.

Prof.

and HOW, Department of Mechani-

cal Egineering, ISMU), Ms. Inderpreet Arora (M.Tech Student, IIT Bombay), Ms. Sandhya Sourirajan (B.E. student, Coimbtore Institute of Technology), all teachers and my friends. I am thankful to the entire work force at ERTS lab, CDEEP and Chemical Department of IIT Bombay for the part they played in making the period of work a joyful experience. Last but not the least, I thank to Indian Institute of Science Banglore (the Organizing Institute) and Department of Science and Technology, Government of India (the sponsoring agency) for their support and encouragement.

3

Abbreviations and meaning of symbols ODE: ordinary dierential equation IVP: initial value problem BVP: boundary value problem

rk: Runge-Kutta 4th order method rkf: Fehlberg's Runge-Kutta order 4 and 5 method (RKF45) LT: Laplace Transform ILT: Inverse Laplace Transforms //**S: Shows the start of the part of program executed online in Scilab window. //**E: Shows the end of the part of program executed online in Scilab window. //##S: Shows start of the writer's contribution to the program or in the matter taken from references or complete work of writer. //##E: Shows end of the writer's contribution to the program or in the matter taken from references or complete work of writer.

Machine congurations Simulation has done in two machines of following congurations 1. Notepad Microsoft Windows Vista Home Premium, Intel(R) Core(TM)2 Duo CPU 1.67GHz 32 bit operating system, 2GB RAM, 160GB HDD. 2. Lab PC Microsoft Windows XP Professional version 2002 SP2, Intel(R) Pentium(R) D CPU 3.40GHz, 1GB RAM, 160GB HDD.

Disclaimer 1. All the programs written in the report are tested in Scilab and produced gure are attached in the report, so in any case not proper working of program is might be due to printing mistake. 2. For purpose of contracting space in the report; position of some plots/gures have changed, however they are correct. 3. In some places in programs '...' this notation has used at last of the lines. Notation is for continuation purpose, in running program this may replace by continuing next line.

4

Contents 1 Scilab

8

1.1

SCILAB HISTORY . . . . . . . . . . . . . . . . . . .

8

1.2

HOW TO OBTAIN ?

8

1.3

INSTALLATION & REQUIREMENTS

. . . . . . . . . . . . . . . . . . . . . . . .

2 Method involve in solution to ODE 2.1

9

Analytical Solutions of Ordinary Dierential Equations (ODEs) 2.1.1

. . . . . . . . . . . . . . . . . . . . . .

9

Introduction to Dierent kind of ODEs and their analytical Solution . . . . . . . . . . . .

2.2

8

9

Numerical methods for solution of ODEs available in Scilab 2.2.1

2.2.2

. . . . . . . . . . . . . . . . . . . . . . . . . .

12

Euler's method (Non-Sti ): Solution of ODEs of non-sti type by the basic denition of dierentiation and by fourier transform expansion. . .

12

Scilab ODE Solvers

12

. . . . . . . . . . . . . .

3 Some Useful Commands regarding to ODE solutions: 13 4 Stability of a solution and Sti/Non-sti problems

13

4.1

Stability . . . . . . . . . . . . . . . . . . . . . . . . .

13

4.2

Sti and Non-sti problems . . . . . . . . . . . . . .

14

5 Scilab Programs for solving ODEs and their examples 16 5.1

Solution to Simple ODE explicit in terms of independent variable

. . . . . . . . . . . . . . . . . . . . . .

16

5.2

Solution by Euler method

. . . . . . . . . . . . . . .

18

5.3

A simple ODE of type dy/dt = f(t,y) . . . . . . . . .

20

5.4

5.5 5.6

Sti type problem and its solution by other methods comparing with sti method . . . . . . . . . . . . . .

21

Solution to Linear rst order linear system of ODEs

25

Solution to the system of Linear ODEs with help of Laplace Transform

5.7

. . . . . . . . . . . . . . . . . . .

Solution to a Second Order ODE with constant coefcients . . . . . . . . . . . . . . . . . . . . . . . . . .

5.8

26

28

Solution to Higher order non-homogeneous ODEs with constant coecients

. . . . . . . . . . . . . . . . . .

5

33

6 Miscellaneous Problems 6.1

36

Predator-Prey model . . . . . . . . . . . . . . . . . . 6.1.1 6.1.2

36

Predator-Prey problem solution by Euler's method 38 Predator-Prey problem solution by Scilab ODE Solver

. . . . . . . . . . . . . . . . . . . . . .

42

6.2

Lorenz ow and Lorenz equation

. . . . . . . . . . .

44

6.3

Curve tting with Scilab for Data analysis purpose .

46

6.4

Solution to ODE where analytical is solution is dicult 51 6.4.1

6.5

Solution of Example 11 by Scialb ode solver .

Modeling and Simulation of a Pressure wave generator

51 52

7 Demerits and diculties in Scilab and suggestion for improvement 55

6

Abstract Techniques and methods for obtaining solutions to dierent kind of Ordinary Dierential Equations is investigated in Scilab.

The

approach is based on solving dierent kind of Ordinary Dierential Equations with dierent method some which are user dened for example Euler's method and other which are ready made in Scilab for example Runge-Kutta, Fehlberg's runge-Kutta, AdamsBashforth and Sti (belonging to sti category problems). Emphasis is placed on mathematical justication of the approach.

Time

required to complete a task and step size for desired accuracy of solution are the main concern and basis of comparison between methods. On the basis of this approach problem and diculties in Scilab are observed and suggestions are made for their remedies.

Preamble In mathematics, an ordinary dierential equation (or ODE) is a relation that contains functions of only one independent variable, and one or more of its derivatives with respect to that variable. Sir Isaac Newton has rst developed the concept of ordinary Dierential Equations for a series solution in the eld of Astronomy and from that time ODE has wide applications in Industry, scientic work, Economics, Ecology, Geology and many other elds including biological sciences. ODE always gives the relation among independent variables, dependent variable and their derivatives which are dicult to analyze. A better analysis of these equations requires their solution in terms of independent and dependent variables. Various techniques are available for the solution of these ODEs in dierent software packages like Mathematica, Maple, Matlab, Scilab etc. This report mainly deals with the solution of ODEs in Scilab by dierent methods like Adams-Bashforth, Runge-kutta , Fehlberg's RungeKutta , Discrete time simulation and nally Sti type of equations along with examples, comparison of their solutions (error involved, smoothness of curve).

7

Aim Aim of this paper is to examine the apparent trend in Simulation of Ordinary Dierential Equation using Scilab by various methods, their applicability and nally analyzing the method with optimize solution for a various category of problems.

1 Scilab 1.1 SCILAB HISTORY It is an open source software package developed at INRIA (France), for system control and signal processing applications.

It also features a wide variety of

tools for various Engineering and Mathematical applications. It was introduced as an Open source alternative to MATLAB. It is also a vector and matrix based program.

Scilab introduced Scicos equivalent to Simulink.

It has constantly

undergone vital changes ever since its inception in 1994.

1.2 HOW TO OBTAIN ? It can be freely downloaded at the home page(http://www.scilab.org/). The site oers the latest version for all the prominent Operating systems. The Windows version is a 83.6 MB zip le ( Scilab 5.1.1). The user is free to modify the source which is also readily available at the same link.

The package also includes a

simulator called SCICOS which is an open source alternative to SIMULINK of MatLab. The minimal package does not include subject specic toolboxes and are to be downloaded separately.

1.3 INSTALLATION & REQUIREMENTS The Scilab installation le has to be unzipped into a folder from which it can be easily installed into the desired location by double clicking the installation le. In case of Linux, the setup le in binary version is a 15.6 MB .tar le which has to be extracted into a directory. Open the Linux terminal and then move to the directory into which the SCILAB .tar le has been extracted.

Typing

the  make command from the SCILAB directory produces an executable le, which can be invoked by typing  scilab .The standard package occupies only 130 MB of disk space upon installation.

8

2 Method involve in solution to ODE 2.1 Analytical Solutions of Ordinary Dierential Equations (ODEs) (Preparation done by attending Continuing Education ProgramCEP02 on Ordinary Dierential Equations at Indian Institute of Technology Bombay held from May 11,09 to May 22,09.)

2.1.1 Introduction to Dierent kind of ODEs and their analytical Solution 1. Simple linear ordinary dierential equation

dy = ax + bxy dt dx = cy + dxy dt 2. Clairut's ODE

y=x

dy dy + f( ) dx dx

3. Singular and General solution, General solution as Envelope of Singular solution. 4. Curl, Divergence, Laplacian of a vector, Convex and Concave domain. 5. Exact or Total Dierential Equation

M dx + N dy = 0 6. Necessary and Sucient condition for Exact or Total Dierential Equation.

∂N ∂M = ∂y ∂x ∂f (x, y) ∂x ∂f (x, y) N= ∂y

M=

7. Along with vector eld

M ˆi + N ˆj

must be conservative in nature.

8. Linear rst order ODE and its solution by operation with Integrating factor.

dy + P (x)y = Q(x) dx

9

9. Bernoulli's Equation and its solution.

dy + P (x)y = Q(x)y n dx 10. Orthogonal Trajectories. 11. Solution of second order Non-homogeneous dierential Equation by method of variation of parameters. 12. Higher order dierential Equation with variable coecients.

Dn y +

n−1 X

aj (x)Dj y = 0

j=0 Calculation of Wronskian. 13. Homogeneous Linear Dierential Equation with constant coecients

Dn y +

n−1 X

aj Dj y = 0

j=0 In cases of Real distinct, Real repeative and complex roots. 14. Non-homogeneous Linear Dierential Equation with constant coecients

Dn y +

n−1 X

aj Dj y = R(x)

j=0 Solution of such problems require Annihilators, which are

T able1 R(x) xk eax xk eax sinbx xk eax cosbx

Annihilator (D − a)k+1 {(D − a)2 + b2 }k+1 {(D − a)2 + b2 }k+1

15. Improper Integrals (1st and 2nd kind). 16. Riemann Integrals and its properties, theorems, tests. 17. Functional function.

f (x + 1) = x(x) and

f (1) = 1

10

18. Gamma function and its properties.



Z

exp(−t)ta−1 dt

Γ(a) = 0 19. Beta function and its properties.

1

Z

xa−1 (1 − x)b−1 dx

β(a, b) = 0 20. Euler's formula

Γ(a)Γ(a − 1) =

π sin(πa)

21. Exponential type function, Necessary and Sucient condition for Laplace Transform. 22. Laplace Transform and its properties.



Z F (s) = L{f (t)} =

exp(−st)f (t)dt 0

23. Shift theorem, Uniqueness Theorem for Laplace Transform. 24. Inverse Laplace Transform , properties, theorems.

L−1 {F (s)} = f (t) 25. Solution of ODE, system of linear dierential equation by Laplace Transform method. 26. Application of Laplace Transform in calculating Integrals. 27. Convolution theorem.

H(s) = F (s) ∗ G(s) then

Z h(t) = (f ∗ g)(t) =

1

f (τ )g(τ − t)dτ 0

28. Homogeneous Linear Dierential Equation with constant coecients and their dierent combinations as:

3

Legendre Equation

3

Tchebychev Equation

3

Airy Equation

3

Hermits Equation

3

Laguerre Equation

3

Bessel's Equation

3

Hypergeometric Equation

3

Jacobi's Equation.

11

2.2 Numerical methods for solution of ODEs available in Scilab 2.2.1 Euler's method (Non-Sti): Solution of ODEs of non-sti type

by the basic denition of dierentiation and by fourier transform expansion.

2.2.2 Scilab ODE Solvers Although there are other ode solvers like ode, dassl, dassrt, odedc are available in Scilab but we will consider only 'ode' solver here. General call syntax of Scilab ODE solver 'ode'

y = ode(0 type0 , y0 , t0 , t, f (t, y)) where y is a function of t and dy/dt = f(t,y)

y0

and

t0 are

initial condition as

y(to ) = y0

t is the range for the solution of ODE type is a character string given in a single quote and direct machine for any specic type of numerical solution method. This methods and corresponding string are as follows;

Default Lsoda method uses non-sti procedure initially and then uses Sti backward dierentiation method , if required.

3

'(empty input)':

3

'adams':

3

'rk':

3

'rkf ':

3

'x':

3

'root':

3

'discrete':

3

'sti ':

Solves Non-sti problems by Adams-Bashforth method.

Uses Runge-Kutta 4th order method to solve non-sti problem.

For non-sti and mildly sti problems, uses Fehlberg's Runge-Kutta pair of order 4 and 5 (RKF45) method, should not use when high accuracy is desired. Similar to Fehlberg's Runge Kutta method, with easier user interface (for non-sti problems). Use default lsoda method and gives roots of the solution (for non-sti problems). For Discrete time simulation, can solve ODE at discrete points (for non-sti problems). Only method which can deal with to deal with sti problems.

Note: However, sti ode solver is particularly for sti problems but other can also solve them under depression of intervals for the whole range of t or for the segment where solution is unstable, same is explained in an example later also.

12

3 Some Useful Commands regarding to ODE solutions: 3

ODE solvers: ode, dassl, dassrt, odedc

3

Polynomials: poly, roots, coe, horner, clean, freq

3

Linear systems: syslin

3

Programming: function, de, argn, for, if, end, while, select,warning, error, break, return

3

Comparison symbols: ==, >=, >, =, & (and),| (or)

3

Execution of a le: exec

3

Character strings: string, part, evstr, execstr

3

Graphics: plot, xset, driver, plot2d, xgrid, locate, plot3d, Graphics

3

Interconnected dynamic systems: scicos

Furthermore, Scilab provides several facilities for nonlinear calculations. Numerical simulation of systems of dierential equations is made by the ode primitive. Many solvers are available, mostly from odepack, for solving sti or non-sti systems. Implicit systems can be solved by dassl. There is a number of optional arguments available for solving ode's (tolerance parameters, jacobian, order of approximation, time steps etc). For ode solvers, these parameters are set by the global variable %ODEOPTIONS. Minimizing non linear functions is done the optim function. Several algorithms (including non dierentiable optimization) are available. Codes are from INRIA's modulopt library. Enter help optim for more a more detailed description. (Reference: http://www.scilab.org/doc/intro/node51.html)

4 Stability of a solution and Sti/Non-sti problems 4.1 Stability Let us take dy/dt = ky, kC solution of this is in the form of y = exp(kt) , this approaches y→

0

as t→



when Re(k)<0. If numerical solution exhibit this

nature then called A-stable solution and don't have stability problems. While solving above equation by Runge-Kutta method we get yn or yn

n

= {ϕ(h, k)} y0 ,

for a stable solution

13

|ϕ(h, k)|

must be <1.

= ϕ(h, k)yn

4.2 Sti and Non-sti problems What's a sti equation? No, precise denition exist. Operationally

dY = f (x, y) dx Y (x0 ) = Y0 is 'Sti ' if its solution by some methods requires (perhaps in a portion of an interval) a signicant depression of the step size in order to avoid instabilities. Take an example of

dy = −15y dt solution of this by Euler's method shows followoing plots for dierent values of steps of t (i.e. h). As h goes smaller solution tends to stable condition.

Although a sti problem can be solve by general methods of numerical solution of ODE with a precision of step increment. Problems other than sti therefore those not show instability while solve by general methods of Solution of ODE are Non-sti problems, also these don't requires more precision on step increment. Another example of sti equation is following linear combination of dierential equations as:

14

and the ODE is

dY = AY dt

solution of this ODE for two values of h; h1 h2

−2

= 2.73972 × 10

= 2.702703 × 10−2 and

is

(References: http://en.wikipedia.org/wiki/stiff_equation,http:// www.physics.arizona.edu/~restrepo/475B/Notes/source/node16.html)

15

5 Scilab Programs for solving ODEs and their examples 5.1 Solution to Simple ODE explicit in terms of independent variable dy = f (x) dx Above can be easily solve by basic integration method like Summation of elemental rectangular area, Trapezoidal rule, Simpon's 1/3 rd etc. Here an examples present to show integration by Summation of elemental rectangular area. Here area between xi and xi+1 can be calculated in three ways a. h*f(xi )

f (xi )+f (xi+1 ) ) 2 c. h*f(xi+1 ) b. h*(

Above three are called Lower sum, Middle sum and Upper sum respectively. Program 1: //##S function [I]=integration(stype,a,b,n,f ) //stype is the way of integration may be 'U','M' & 'L' for upper, middle & lower sum //a,b is the range of x for integration, n=no. of intervals, f be the function dened online h=(b-a)/n; xset('window',1) if (stype<>'L')&(stype<>'M')&(stype<>'U') then error('stype must be L,M,U as string '); abort; end if stype=='L' then x=a & j=1 elseif stype=='M' then x=a+(h/2) & j=2 else x=a+h & j=3 end A=zeros(1,n) X=zeros(1,n) Y=zeros(1,n) for i=1:n X(1,i)=x y=feval(x,f ) Y(1,i)=y

16

A(1,i)=h*y x=x+h end subplot(2,2,1) plot(X,Y) title('x,y function plot') I=sum(A) B=zeros(1,n); for i=1:n B(1,i)=sum(A(1,1:i)) end subplot(2,2,2) plot(X,B) title('Integrated curve for '+stype) endfunction //##E The function name above is given as integration and calling syntax is integration(stype,a,b,n,f ), here stype is the type of integration (Lower, Middle or Upper sum), a and b are limit of integration, n is the no. of intervals for acccuracy of result. Example 1: Let

dy = x2 + x + 3 dx

//**S //##S de('[z]=f(x)','z=x^2+x+3'); getf('integration.sci'); integration('L',0,5,10,f ) ans = 68.41875 integration('M',0,5,10,f ) ans = 69.165625 integration('M',0,5,10,f ) ans = 69.91875 //**E //##E Three answer shown above are for Lower, Middle and Upper sum respetively, let them call L-ans, M-ans and U-ans. Note:

17

L − ans < M − ans < U − ans and from the exact solution we have answer as 69.166667, which is closest to the M-ans. Why? This can be explained by gure. Fig 1.4 shows the plot of integrated curve of given problem, there are three curves for each sum and so make the curve line thick.

Fig 1 Note: However, above program is correct to get answer and individual plots but given Fig 1 is not produce by above program. Some more programming is required for the same.

5.2 Solution by Euler method as we know of the ODE is

dy = g(t, y) dt

then

g(y + h) = g(y) + hg(t, y) 18

or

yn+1 = yn + hg(tn , yn ) Following above program for solution by Euler's method can be written as: Program 2: function [t,y]=Euler(t0,y0,tn,h,g) //Euler 1st order method for solving ODE //dy/dx=g(t,y) //t0 and tn are the range of t and h is the interval bet ymaxallowed=1e+100; //Actually it is abs(ymaxallowed) t=t0:h:tn; y=zeros(t); n=length(y); y(1)=y0; for j=1:n-1 y(j+1)=y(j)+h*g(t(j),y(j)); if abs(y(j+1))>ymaxallowed then disp('Euler-warning: underow or overow') n=j; t=t(1,1:n); y=y(1,1:n); break; end; end; endfunction;

(Reference: numerical and statistical methods with SCILAB for science and engineering, vol. 1, gilberto e. urroz) Example 2: //Let the function is

dy dt

= f (t, y) =

1 2(1+y) with initial condition y(0)=-1

,from t:2 to 10 //**S //##S t0=1; y0=0; de('[z]=g(t,y)','z=1./(2*(1+y))'); getf('Euler.sci') //calling earlier dened Euler program Euler(t0,y0,tn,h,g); plot(t,y) xtitle('Fig 2-Solution plot to Example 1 by Euler method','t','y') //##E //##S

19

Fig 2

5.3 A simple ODE of type dy/dt = f(t,y) Example 3: Let us solve the problem of Example 2 again //**S //##S >t0=1; y0=0; //initial condition >de('[z]=f(t,y)','z=1./(2*(1+y))')//dening function >t=1:1:10; >y=ode(y0,t0,t,f ); >plot(t,y) >xtitle('Fig 1-Solution plot to Example 1','t','y') //**E //##S It can be observe that the output plots between dependent and independent variable are similar for both the method therefore by Euler's method and Scilab ODE solver. It is dicult to distinguish them with such simple problems. Their performance may be distinguish in case of complicated problems or in system of ODEs with large no. of iterations.

20

Fig 3

5.4 Sti type problem and its solution by other methods comparing with sti method What is Sti type of problem? Example 4: Let us take an example dy/dx=-1000(y-x)+2001 (whose exact solution is y1 (x)=-exp(-1000x) and y2 (x)=x+2 ) We will plot a particular solution y=y1(x)+y2(x)=-exp(-1000x)+x+2 of given ODE Note: In particular solution rst part goes zero rapidly while second part produces a straight line Program 4: //##S de('[y]=f(x)','y=-exp(-1000.*x)+x+2') //PHASE 1: Plotting exact solution for small values of x x1=linspace(0,0.01,100); y1=feval(x1,f ); xset('window',0) clf subplot(4,1,1) plot(x1,y1)

21

xtitle('Fig 4.1:Plot for exact solution for small values of x','x','y') //PHASE 2:Plotting exact solution for small values of x x2=linspace(0,1,100); y2=feval(x2,f ); subplot(4,1,2) plot(x2,y2) xtitle('Fig 4.2:Plot for exact solution for comparatively large values of x','x','y') //PHASE 3:Now we will see what happen when we plot the curves for dierent numerical solutions de('[Dy]=g(x,y)','Dy=-1000.*(y-x)+2001') //Call 'Euler' program for Numerical solution by Euler method getf('Euler.sci') //Plot the curve for dierent increment as 0.0005,0.0010,0.0020,0.0025 //For initial condition y(0)=1 y0=1; x0=0; h1=0.0005; h2=0.0010; h3=0.0020; h4=0.0025; [x1,y1]=Euler(x0,y0,0.01,h1,g); [x2,y2]=Euler(x0,y0,0.01,h2,g); [x3,y3]=Euler(x0,y0,0.01,h3,g); [x4,y4]=Euler(x0,y0,0.01,h4,g); ymin=min([y1 y2 y3 y4]); ymax=max([y1 y2 y3 y4]); disp(ymin) disp(ymax) //Let us dene the axis properties rect=[0 -4 0.01 6]; subplot(4,1,3) //plots to give discrete curve with dierent signs plot2d(x1,y1,-1,'011',,rect) plot2d(x2,y2,-2,'011',,rect) plot2d(x3,y3,-3,'011',,rect) plot2d(x4,y4,-4,'011',,rect) //plots to give continuous curves plot2d(x1,y1,1,'011','h=0.0005',rect) plot2d(x2,y2,5,'011','h=0.0010',rect) plot2d(x3,y3,3,'011','h=0.0020',rect) plot2d(x4,y4,2,'011','h=0.0025',rect) title('Fig 4.3:Plots by Euler method for dierent value of h (increment)') //Here corresponding curves for dierent colors are as follows //h=0.0005, color=Black //h=0.0010, color=Red //h=0.0020, color=Pale Green //h=0.0025, color=Blue //PHASE 4:Now Let us see what happen when we plot the above curves with sti ODE solver of Scilab subplot(4,1,4)

22

rect=[0 1 0.01 2.2] x1=x0:h1:0.01; x2=x0:h2:0.01; x3=x0:h3:0.01; x4=x0:h4:0.01; y1=ode(y0,x0,x1,g); y2=ode(y0,x0,x2,g); y3=ode(y0,x0,x3,g); y4=ode(y0,x0,x4,g); plot2d(x1,y1,1,'011','h=0.0005',rect) plot2d(x2,y2,5,'011','h=0.0010',rect) plot2d(x3,y3,3,'011','h=0.0020',rect) plot2d(x4,y4,2,'011','h=0.0025',rect) title('Fig 4.4:Plot for dierent h (x-increment) by sti solver of Scilab') //From the plot it is obvious that sti ODE solver of Scilab stabilizes solution for x-increment=0.0005,0.0010,0.0020,0.0025 //##E Answer: - 3.0525, 5.3825 Explanation: However, program itself explains everything but still detailed explanation is as follows: Fig 4.1 is the plot of exact solution

{y(x) = −exp(−1000x)andy(x) = x + 2}

for the smaller values of x. Fig 4.2 gives the same solution plot for comparatively larger values of x. It is obvious from the plots that solution faces sudden changes between 0.0 and 0.01. Fig 4.3 gives the solution of ODE by Euler's Method at dierent values of x-increment (h), for those curves are shown with dierent colors. Answer given after execution of program is the minimum and maximum value of y respectively, comes in the solution of ODE by Euler's method (just to check degree of instability). Color and their corresponding (h) curves are as follows: Table 2 Value of regular x-increment

Color of curve

h=0.0005

Black

h=0.0010

Red

h=0.0020

Pale Green

h=0.0025

Blue

Apparently as value of interval increases (h) instability of curve increases.

23

Fig 4.4 comes after solution of given ODE by SCILAB Sti ODE solver at all previously dened values of h (colors of curve are corresponds to values of h as similar to Fig 4.3). It is apparent from the plots that sti solver stabilize the solution but in better way for higher value of h.

Fig 4

(Reference: numerical and statistical method with SCILAB for science and engineering, Volume 1, gilaberto e. urroz) 24

5.5 Solution to Linear rst order linear system of ODEs This is example gives the idea behind how to solve system of ODE with a single ode solver Example 5: Let us take system of ODEs as

dy1 = y2 + x dx and

dy2 = −y1 + y2 dx

with initial condition as y1 (0)

= 1;y2 (0) = 2;between

0 to 2 with regular increment of 0.1

Program 5: //**S de('[w]=f(x,y)',['f1=y(2)+x';'f2=-y(1)+y(2)';'w=[f1;f2]']) x0=0; Dx=0.1; xn=2; y0=[1;2]; x=[x0:Dx:xn]; y = ode(y0,x0,x,f ); plot2d([x',x'],[y(1,:)',y(2,:)'],[1,-1],'111','y1@y2',[0 -3 2 4]) xtitle('ode solution to a system of ODEs in example 3','x','f(x)') pause //**E In the Fig 5 continuous plot gives the solution in discrete pattern shows solution

y2 .

25

y1

and the plot which is shown

Fig 5

(Reference: Youngstown state University,http://www.eng.ysu. edu/~jalam/engr6924s07/sessions/session27/session27.pdf)

5.6 Solution to the system of Linear ODEs with help of Laplace Transform Let a system is as

dY = A ∗ Y + g(t) dt where A is a n×n matrix if n no. of ODEs are there. Taking Laplace Transform at both side

L{

dY } = A ∗ L{Y } + L{g(t)} dt

sL{Y } − Y0 = A ∗ L{Y } + L{g(t)} (sI − A)L{Y } = Y0 + L{g(t)} L{Y } = (sI − A)−1 Y0 + (sI − A)−1 L{g(t)} 26

Taking Inverse Laplace Transform at both side

Y = L−1{ (sI − A)−1 }Y0 + L−1 {(sI − A)−1 L{g(t)}} Program 6: Following program has been made in MatLab. //##S function [Y]=LT(A,x0) [m,n]=size(A); if m>n error('matrix must of square type') elseif m
dy1 = 4y1 + 2y2 dt dy2 = −2y2 dt for initial conditions as y1 (0)

=0

and y2 (0)

This problem can be re-written as

dY =A∗Y dt where

A=

4 0

//**S //##S A=[4 2;0 -2]; x0=[0 1]; LT(A,x0) //##E //**E

27

2 −2

=1

and 0
Answers: ParticularSolutions = y1=2/3*exp(t)*sinh(3*t) y2=exp(-2*t) GeneralSolutions = y1=exp(4*t)*C1+2/3*exp(t)*sinh(3*t)*C2 y2=exp(-2*t)*C2

5.7 Solution to a Second Order ODE with constant coefcients This type of ODE is discussed here, because it comes frequently in dierent disciplines of Engineering especially in Electrical, Mechanical and Civil Engineering. Let us take an example of spring-mass system as shown

Let m be the mass, k be the spring constant and b be the damping constant which appear due to frictional losses (may be heat losses in case of Electrical LCR circuits). And forces corresponding to above are

F = −kx F = −bv 28

where x and v are instant position and velocity of mass respectively. From Newton's

2nd law

−kx − bv = m

−kx − b

d2 x dt2

dx d2 x =m 2 dt dt

d2 x k b dx + x+ =0 2 dt m m dt D2 +

b k D+ =0 m m

which is a quadratic equation whose roots may be real distinct, real repeated and complex roots. A particular type of roots gives the relation among m, b and k or inequalities among them which is important because a particular type of root decides the characteristic of system and its functioning, this can be seen as follows: Program to solve above ODE Program 7: //##S function []=dampedoscillation(m,B,k,t0,x0,v0,t1,t2) //m=mass, b=damping constant, k=spring constant, t0=initial time... //at which xo,vo are given, //t=time range for plot, //F1=-kx, F2=-bv //Di. Equation: (D^2+(b/m)D+(k/m))x=0 t=linspace(t1,t2,1000); D=poly([(k/m) (B/m) 1],'D','coe ') R=roots(D) disp('roots') disp(R) a=R(1,1) b=R(2,1) //PHASE 1: For distinct real roots if imag(a)==0 & a<>b then M=[exp(a*t0) exp(b*t0);a*exp(a*t0) b*exp(b*t0)]\[x0 v0]' //A & B are coecients in the solution x=Aexp(R1*t)+Bexp(R2*t) A=M(1,1)

29

B=M(2,1) x=(A*exp(a*t))+(B*exp(b*t)) v=(a*A*exp(a*t))+(b*B*exp(b*t)) a=(a*a*A*exp(a*t))+(b*b*B*exp(b*t)) clf subplot(2,2,1) plot(t,x) title('x-t curve') subplot(2,2,2) plot(t,v) title('v-t curve') subplot(2,2,3) plot(t,a) title('a-t curve') subplot(2,2,4) plot(t,x,t,v,t,a) title('x-t, v-t & a-t curve at one place') //PHASE 2: For repeated real roots elseif a==b then M=(exp(a*t0))*[1 t0;a (a*t0)+1]\[x0 v0]' //A & B are coecients ... //in the solution x=Aexp(R1*t)+Bexp(R2*t) A=M(1,1) B=M(2,1) x=(A*exp(a*t))+(B*t.*exp(a*t)) v=(a*A*exp(a*t))+(B*(a*t+1).*exp(b*t)) a=(a*a*A*exp(a*t))+ ...(a*B*(a*t+2).*exp(b*t)) clf subplot(2,2,1) plot(t,x) title('x-t curve') subplot(2,2,2) plot(t,v) title('v-t curve') subplot(2,2,3) plot(t,a) title('a-t curve') subplot(2,2,4) plot(t,x,t,v,t,a) title('x-t, v-t & a-t curve at one place') //PHASE 3: For complex roots else b=imag(a) b=abs(b) a=real(a) c=(a*cos(b*t0))-(b*sin(b*t0)) d=(a*sin(b*t0))+(b*cos(b*t0))

30

M=(((exp(a*t0)))*[cos(b*t0) sin(b*t0);c d])\[x0; v0] //A & B are coecients in the solution x=Aexp(R1*t)+Bexp(R2*t) A=M(1,1) B=M(2,1) x=((A*cos(b*t))+(B*sin(b*t))).*exp(a*t) v=(exp(a*t)).*(((a*A+b*B)*cos(b*t))+((a*B-b*A)*sin(b*t))) a=(exp(a*t)).*((((a*a*A+b*a*B)*cos(b*t))+((a*a*B-b*a*A)*sin(b*a*t)))+... ...(((-b*a*A-b*b*B)*sin(b*t))+((a*b*B-b*b*A)*cos(b*a*t)))) // Note that in all above three equation 'a' i.e. the real part of the roots has the role for damping similarly 'b' i.e. the imaginary part of roots... //has the role for frequency or say angular frequency=b clf subplot(2,2,1) plot(t,x) title('x-t curve') subplot(2,2,2) plot(t,v) title('v-t curve') subplot(2,2,3) plot(t,a) title('a-t curve') subplot(2,2,4) plot(t,x,t,v,t,a) title('x-t, v-t & a-t curve at one place') end xset('window',5) clf subplot(3,1,1) plot(x,v) title('v(vertical axis)-x(horizontal axis)') subplot(3,1,2) plot(v,a) title('a(vertical axis)-v(horizontal axis)') subplot(3,1,3) plot(x,a) title('a(vertical axis)-x(horizontal axis)') endfunction; //##E Now we will see characteristics of system by their response plots for three cases

3

Real distinct roots

3

Real repeated roots

3

Complex roots

31

Fig 6 Explanation: To plot above gures following of m, b and k have been taken Table 3 Root type

m

b

k

Roots -0.2928,-1.7071

Real distinct

2

4

1

Real repeated

2

4

2

-1,-1

Complex

2

4

50

-1±4.8989795 i

32

Assume Fig 6 a matrix of 4×3 then the rst row plots the position vs time (x-t), velocity vs time (v-t) and acceleration vs time (a-t). Second, third and forth rows gives velocity vs position (v-x), acceleration vs velocity (a-v) and acceleration vs position (a-x) plots respectively. Conclusions from the plot:

3

If damping constant is not zero then all motion decays.

3

In case of both real root cases motion is only decaying not oscillatory.

3

To make motion oscillatory roots must be complex which concludes

b2 < 4mk this gives an inequlity between damping constant, mass and spring constant.

3

In the complex root case when damping constant b6=

0, oscillatory motion

goes on reducing which is obvious from all last plots of rst, second and third rows.

3

In case of oscillatory motion acceleration is always proportionate to position of massive body.

5.8 Solution to Higher order non-homogeneous ODEs with constant coecients Dn y +

n−1 X

aj Dj y = R(x)

j=0 or

Dn y = −

n−1 X

aj Dj y + R(x)

j=0 Example 7: Let the problem is

d4 y d3 y d2 y dy x2 + 3 − 2 + 5 + y = dx4 dx3 dx2 dx 2 or

d4 y d3 y d2 y dy x2 = −3 3 + 2 2 − 5 −y+ 4 dx dx dx dx 2

for initial conditions at x=0

D3 y = −1, D2 y = 0, Dy = −1,y = 1

33

This problem can be re-written as

  −3 u3 (x)   d  u (x)  2 = 1    0 u (x) dx 1 0 y(x) 

where

ui (x)are

2 −5 0 0 1 0 0 1

  u3 (x) −1  u2 (x) 0  ∗ 0   u1 (x) y(x) 0

 0.5x2   0   +   0  0 



intermediate function of x.

for this initial condition matrix can be written as



 −1  0   u(0) =   −1  1 Now above problem can be easily solve by Scilab ode solver for 0
34

xgrid(5) xtitle('solution to the 4th order non-homogeneous ODE... (constant coecient), Abscissa-time',,'(D^3)y') yprim4=-3*v(1,:)+2*v(2,:)-5*v(3,:)-v(4,:)+(t.^2)/2; subplot(5,1,5) plot(t,yprim4) xgrid(5) xtitle('solution to the 4th order non-homogeneous ODE... (constant coecient), Abscissa-time',,'(D^4)y') pause

solution to the 4th order non−homogeneous ODE (constant coefficient), Abscissa−time

y

120 100 80 60 40 20 0 −20

0

1 2 3 4 5 6 7 8 9 solution to the 4th order non−homogeneous ODE (constant coefficient), Abscissa−time

10

0

1 2 3 4 5 6 7 8 9 solution to the 4th order non−homogeneous ODE (constant coefficient), Abscissa−time

10

0

1 2 3 4 5 6 7 8 9 solution to the 4th order non−homogeneous ODE (constant coefficient), Abscissa−time

10

0

1 2 3 4 5 6 7 8 9 solution to the 4th order non−homogeneous ODE (constant coefficient), Abscissa−time

10

Dy

100 80 60 40 20 0 −20

(D^2)y

60 40 20 0 −20 −40 −60 −80

(D^3)y

40 20 0 −20 −40 −60 −80 −100 −120 −140 −160

(D^4)y

40 20 0 −20 −40 −60 −80 −100 −120

0

1

2

3

4

5

6

7

8

9

10

Fig 7

(Reference: numerical and statistical method with SCILAB for science and engineering, Volume 1, gilaberto e. urroz)

35

6 Miscellaneous Problems 6.1 Predator-Prey model Model describes the interaction of two species in an eco-system when they are isolated from others. Here we will discuss growth rate of Fish and Shark when they don't interact with other species (assumed). In this situation rate changes will follow relations given below Growth rate of sh=Rate at which sh born-Rate at which sh are eaten by sharks Growth rate of shark=Rate at which sh turned into sharks-Rate at which shark die without sh Now Rate at which sh born∝F Rate at which sh are eaten by sharks∝FS Rate at which sh turned into sharks∝FS Rate at which shark die without sh∝S where F and S are instant population of Fish and Shark respectively. So, we have

dF = αF − βF S dt

(1)

dS = βF S − γS dt

(2)

with the initial conditions

F (0) = F0 S(0) = S0 where F0 and S0 are initial population at time t=0 (relatively).

α :growth

rate of sh in absence of shark (1/years)

γ :death

rate of sharks in the absence of their prey i.e. sh (1/years)

β :death

rate per encounter of sh with sharks (1/sharks/years)

36

 :eciency

of turning predated sh into shark (shark/sh)

From above equations, it is apparent that 1. In absence of shark

dF = αF dt F = F0 eαt

(3)

Thus population of sh increases exponentially. 2. In absence of sh

dS = −γS dt S = S0 e−γt

Thus population of sh decreases exponentially.

Non-dimensionalization: Non-dimensionalization of equations is a process to make all terms dimensionless. Benets of this are described later If we nondimensionalize according to

F∗ =

F F0

S∗ =

S S0

t∗ = αt putting this in equations 1 and 2, we get

dF ∗ βS0 ∗ ∗ = F∗ − F S ∗ dt α dS ∗ βF0 ∗ ∗ γ ∗ F S − S = dt∗ α α if we write

a=

βS0 α

b=

βF0 α 37

(4)

c=

γ α

Then equations can be re-writtem as

dF ∗ = F ∗ − aF ∗ S ∗ dt∗

(5)

dS ∗ = bF ∗ S ∗ − cS ∗ dt∗

(6)

with initial conditions

F ∗ (0) = 1 S ∗ (0) = 1 So, benets of non-dimensionalization Equation becomes easy by reducing no. of undetermined constants. Initial condition are not undetermined for solution of this ODE.

6.1.1 Predator-Prey problem solution by Euler's method Program 9: Program for the Solution of the above problem. //##S function []=PP(Alpha,Beta,Gamma,Epsilon,t,Dt,F0,S0) //Note: t must be in years (absolute, not by transformation) and... //Dt must be in terms of months a=Beta*S0/Alpha; b=Epsilon*Beta*F0/Alpha; c=Gamma/Alpha; t=t*Alpha; Dt=(1/52)*Dt; //Since Dt is given in terms of weeks, so converting it to years Dt=Alpha*Dt; t=0:Dt:t; //t here is the vector from 0 to t multiplied by Alpha //We have {dF^{*}}/{dt^{*}}=F^{*}-aF^{*}S^{*} and... //{dS^{*}}/{dt^{*}}=bF^{*}S^{*}-cS^{*} as transformed equations //Above equations can be re-written as dy1/dt=y1-a*y1*y2 and... //dy2/dt=b*y1*y2-c*y2 //Now the RHS of both above equations are assumed as f1 and f2 d=size(t); e=d(1,2); f1=zeros(1,e);

38

f2=zeros(1,e); y1=zeros(1,e+1); y2=zeros(1,e+1); y1(1,1)=1; y2(1,1)=1; f1(1,1)=y1(1,1)-a*y1(1,1)*y2(1,1); f2(1,1)=b*y1(1,1)*y2(1,1)-c*y2(1,1); for i=1:e f1(1,i+1)=y1(1,i)-a*y1(1,i)*y2(1,i); f2(1,i+1)=b*y1(1,i)*y2(1,i)-c*y2(1,i); y1(1,i+1)=y1(1,i)+Dt*f1(1,i); y2(1,i+1)=y2(1,i)+Dt*f2(1,i); end y1max=max(y1); y2max=max(y2); y1min=min(y1); y2min=min(y2); Fmax=y1max*F0 Smax=y2max*S0 Fmin=y1min*F0 Smin=y2min*S0 disp(Fmin, Fmax, Smin, Smax) disp(Fmin) disp(Fmax) disp(Smin) disp(Smax) if y1min < y2min then p=-y1min else p=-y2min end if y1max > y2max then q=y1max+0.1*y1max else q=y2max+0.1*y2max end r=(max(t))/Alpha; rect=[0 p r q]; t=t/Alpha; subplot(2,1,1) plot2d(t,y2(1,1:e),5,'111',,rect) plot2d(t,y1(1,1:e),11,'111',,rect) xtitle('Fish and Shark population Vs time','time',... 'Shark(S/S0-Red),Fish(F/F0-Blue)') subplot(2,1,2) plot(y1(1,1:e),y2(1,1:e))

39

xtitle('Shark Vs Fish relative population plot',... 'Fish(F/F0)','Shark(S/S0)') endfunction //##E Example 8: Following are the results come from the solution by Euler's method for given values of

α = 0.7, β = 0.007, γ = 0.5,  = 0.3, t = 50, F0 = 200, S0 = 50 Table 4 Time step(weeks)

Fmin

Fmin

Smin

Smax

1.00e+00

58

654

30

226

5.00e-01

79

543

40

206

2.50e-01

90

501

45

191

1.25e-01

95

483

47

184

6.25e-02

98

474

48

181

3.13e-02

99

470

45

180

1.56e-02

100

468

49

179

7.81e-03

100

467

49

178

Conclusion: 1. As steps goes smaller, answers tends to stabilized values. 2. It is aapprent from the plots that as time steps goes smaller solution tends to cyclic solution with invariable cycle amplitudes. Here, Fish and Shark relative population plot with time and their phase portrait for three values of time steps 1.00e+00, 6.25e-02 and 7.81e-03 are given.

40

Fig 8

However by Euler' method solution of Predator-Prey problems are possible but for accurate result (i.e. cyclic curve Population Vs time, should be of constant amplitudes or Predator Vs Prey population plot should be thin in nature) time step should be suciently less and while doing such simulation, it takes a longer time than usual and sometimes due to limited stack size not gives the desired result.

(Reference: Stanford University, http://fluid.stanford.edu/ ~finger/teaching/numerical_methods_02/tutorials/tutorial2.pdf)

41

6.1.2 Predator-Prey problem solution by Scilab ODE Solver Program 10: //##S function []=PPS(stype,Alpha,Beta,Gamma,Epsilon,t,Dt,F0,S0) //Here stype is the method by which... //solution is going to be taken and it may be 'adams', 'rk', 'rkf ', 'sti ' etc. a=Beta*S0/Alpha; b=Epsilon*Beta*F0/Alpha; c=Gamma/Alpha; t=t*Alpha; Dt=(Dt/52)*Alpha; t=0:Dt:t; y0=[1;1]; t0=0; de('[w]=f(x,y)',['f1=y(1)-a*y(1).*y(2)';'f2=b*y(1).*y(2)-c*y(2)';'w=[f1;f2]']) y=ode('stype',y0,t0,t,f ); t=t/Alpha; subplot(2,1,1) plot2d([t',t'],[y(1,:)',y(2,:)'],[11,5],'111','Fish@Shark',[0 -0.5 50 5]) title('Fish and Shark population Vs time') subplot(2,1,2) plot(y(1,:),y(2,:)) xtitle('Shark Vs Fish relative population plot','Fish(F/F0)','Shark(S/S0)') endfunction //##E Above Predator-Prey problem has been also solved by adams, rk, rkf , sti method and conclusions after them are as follows: 1. All these methods are able to give ne results even for a time step of 1 week. 2. Problem is not of sti type because for a same time step solution by method other than sti method and solution by sti method are similar. Solution is already stabilize in nature so there is no need for sti solution. 3. It has been observed that for larger time steps (for eg: 1 year) plots are non smooth and similar for all method, one of this i.e. by rk method is given in Fig 10. 4. Plots in Fig 9 are given for adams, rk, rkf and sti method with time step of 1 week.

42

Fig 9.1

Fig 9.2

43

Fig 10

6.2 Lorenz ow and Lorenz equation The Lorenz equation was published in 1963 by a meteorologist and mathematician from MIT called Edward N. Lorenz.This is a model for some of the unpredictable behavior which we normally associate with the weather. The Lorenz equation is commonly dened as three coupled ordinary dierential equation like

dx = σ(y − x) dt dy = x(τ − z) − y dt dz = xy − βz dt where

44

where the three parameter

σ, τ , β

are positive and are called the Prandtl

number, the Rayleigh number, and a physical proportion, respectively. It is important to note that the

x, y , z

are not spacial coordinate. The  x is

proportional to the intensity of the convective motion, while

y

is proportional

to the temperature dierence between the ascending and descending currents, similar signs of

x

and

y

denoting that warm uid is rising and cold uid is

descending. The variable

z

is proportional to the distortion of vertical

temperature prole from linearity, a positive value indicating that the strongest gradients occur near the boundaries. Program 11: //##S function []=Lorenz(stype,Sigma,Tou,Beta,x0,y0,z0) //Let for convenience xset('window',4) clf a=Sigma; b=Tou; c=Beta; de('[w]=f(t,u)','w=[a*(u(2)-u(1)); u(1).*(b-u(3))-u(2); u(1).*u(2)-c*u(3)]'); t0=0; t=0:0.01:50; u0=[x0 y0 z0]'; u=ode(stype,u0,t0,t,f ); subplot(3,1,1) plot(t,u(1,:),,'t','x') subplot(3,1,2) plot(t,u(2,:),,'t','y') subplot(3,1,3) plot(t,u(3,:),,'t','z') xset('window',2) clf plot3d3(u(1,:),u(2,:),u(3,:)) title('Phase portrait of x,y,z') endfunction //##E Example 9: Lorenz Equation is famous for its drastic dependency on initial conditions and parameters

σ, τ, β .

We will take here three examples by changing parameters to see their eect. 1.

σ

= 10,

τ

= 28 and

β

= 8/3 (which is the classical example).The initial

condition of the system is (x0 , y0 , z0 )= (3,15,1). 2.

σ

= 5,

τ

= 28 and

β

= 8/3 (which is the classical example).The initial

condition of the system is (x0 , y0 , z0 )= (3,15,1).

45

3.

σ

= 10,

τ

= 100 and

β

= 8/3 (which is the classical example).The initial

condition of the system is (x0 , y0 , z0 )= (3,0,1). Conclusions: 1. For all rk, rkf, adams solution are similar. 2. Fig 11 shows how solution changes drastically with variation in parameter and initial conditions.

Fig 11

(Reference: http://planetmath.org/encyclopedia/lorenzEquation.html)

6.3 Curve tting with Scilab for Data analysis purpose Curve tting is nding a curve which has the best t to a series of data points and possibly other constraints. This section is an introduction to both interpolation (where an exact t to constraints is expected) and regression analysis. Both are sometimes used for extrapolation. Regression analysis allows for an approximate t by minimizing the dierence between the data points and the curve. Technique to t set of data in a curve (or polynomial) is based on minimizing their sum of square of dierence between experimental data and data provided by considered polynomial.

46

Here a program has been made without using Scilab readymade curve tting functions. Program can t up to ve curve at a time and shoes their plot at one place for analysis purpose. Function curvetting given below requires three input these are vector of independent variable, vector of experimental data and a vector which may have size up to ve and contains degree of polynomials on which curve has to be tted. As the result of executing program it gives the coecient of individual terms of a polynomial. Program also gives the error plot for each tted curve and value of sum of squares of dierence between experimental and value provided by plot.

y = a0 + a1 x + ax22 + a3 x3 + a4 x4 ........ Program 12: //##S //let k be the degree of the polynomial //Fundamental equation A*a=B function []=curvetting(x,y,D) xset('window',1); clf //let us rst dene how many elements are there in D p=size(D); q=p(1,2); //##PHASE1: Getting A subplot(6,2,1); plot(x,y) xgrid(#) title('DATA CURVE') m=size(x); //n actually gives the no. of elements in the vector x or y n=m(1,2); for l=1:q // starting of this for loop k=D(1,l) A=zeros(k+1,k+1); for j=1:k+1 for i=1:k+1 A(i,j)=sum(x.^(i+j-2)) end end //##PHSAE2:Getting B B=zeros((k+1),1); for i=1:k+1; B(i,1)=sum(y.*(x.^(i-1))) end //##PHASE3:Getting 'a' the coecient matrix

47

//'a' is the desired coecient matrix shows a0,a1,a2... // in the tting polynomial a=(A\B)'; disp(a) //##PHASE4:Ploting curve by dierent solutions //c & d be the rst & last value of vector x c=x(1,1); d=x(1,n); x1=linspace(c,d,50); //getting the smooth value of x x2=x1'; //getting transpose of x1 //we have made x1 because to get a smooth curve of resultant polynomial x3=zeros(50,k+1); // an assumed matrix for further operation //we going to make x3 a matrix having rst row as... //1,x2.^2,x2.^3.... for rst value of x2, similarly second row... //will be for second value of x2, so in the above manner we... //will get a matrix having its row as 1,x2.^2,x2.^3.... //now if any of the row of this matrix multiply with 'a'... //& suppose then give a vector then sum of that vector... //element will be the corresponding y2(or say y3 here) for x2 x3(:,1)=1; x3(:,2)=x1'; for i=1:k+1; x3(:,i)=x2.^(i-1) end y3=zeros(50,1); for i=1:50 y3(i,1)=sum(a.*x3(i,:)) end //dening y4 as calculated from result polynomial according to 'a' matrix subplot(6,2,(2*(l-1)+3)) plot(x2,y3) xgrid(3) title('CURVE FOR DEGREE='+string(k)); //PHASE5:To get the value of y corresponding to given x for resultant polynomial x4=zeros(n,k+1); for i=1:k+1; x4(:,i)=x'.^(i-1) end y4=zeros(n,1); for i=1:n y4(i,1)=sum(a.*x4(i,:)) end

48

//PHASE6:To get the error curve for corresponding plot,... //i.e. plot between (ye-yr) & x. where ye is the experimental... // value and yr is the real value z=y'-y4; z1=z'.^2; z2=sum(z1); z3=sqrt(z1); b=[k,z2]; subplot(6,2,2*(l+1)) plot(x,z3) xgrid(3) title('SUM OF SQUARE OF DEVIATIONS FOR THIS ERROR CURVE='+string(z2)); ylabel('DEGREE OF CURVE='+string(k)); end // end of rst for loop endfunction; //##E Example 10.1: Let us try to t the known data for a circle of radius 5 centered at origin. Producing data //**S >de('[y]=f(x)','y=sqrt(25-x.^2)') >x=-5:0.1:5; >y=feval(x,f ); //**E Now tting the generated data for polynomial of degrees 1,2,3,4, and 5 by using the program 12. It is apparent from the plots second and third degree curves are more or less representing the circle. The answer came from program are as follows: Degree

Coecients

Sum of square of deviations

1

3.9,-1.6D-17

142.87233

2

5.2,1.0D-16,-0.2

7.0324985

3

5.2,-2.5D-17,-0.2,8.5D-18

7.0324985

4

5.0,2.0D-16,-0.1,-3.0D-18,-0.01

1.8451398

5

5.0,-6.0,-0.1,1.1D-15,-0.01,-4.0D-17

1.8451398

49

Fig 12.1 From the above plot conclusion can be made that while approaching to higher degree plots although sum of square of deviation (sum of least squares) reduces but along with this smoothness of the curve is also reduces, which are opposite to each other for a desired result in industry. So one should go with the optimize result. Example 10.2: Similar thing has been done a Gaussian curve for

σ = 0.9 ,xmean = 5and

thus the Fig 12.2 came out concludes that no tted curve able to represent the Gaussian curve.

50

DATA CURVE

−5

0 5 10 CURVE FOR DEGREE=1

15

20

−5

0 5 10 CURVE FOR DEGREE=2

15

20

−5

0 5 10 CURVE FOR DEGREE=3

15

20

−5

0 5 10 CURVE FOR DEGREE=4

15

20

−5

0 5 10 CURVE FOR DEGREE=5

15

20

15

20

−5

0

5

10

DEGREE OF CURVE=5 DEGREE OF CURVE=4 DEGREE OF CURVE=3 DEGREE OF CURVE=2 DEGREE OF CURVE=1

0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 −10 1.0 0.8 0.6 0.4 0.2 0.0 −0.2 −0.4 −0.6 −0.8 −1.0 −10 0.08 0.06 0.04 0.02 0.00 −0.02 −0.04 −0.06 −10 0.08 0.06 0.04 0.02 0.00 −0.02 −0.04 −0.06 −10 0.03 0.02 0.01 0.00 −0.01 −0.02 −0.03 −0.04 −0.05 −10 0.015 0.010 0.005 0.000 −0.005 −0.010 −0.015 −0.020 −10

SUM OF SQUARE OF DEVIATIONS FOR THIS ERROR CURVE=2.8021607 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 −10 −5 0 5 10 15 20 SUM OF SQUARE OF DEVIATIONS FOR THIS ERROR CURVE=2.3957277 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 −10 −5 0 5 10 15 20 SUM OF SQUARE OF DEVIATIONS FOR THIS ERROR CURVE=2.3957277 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 −10 −5 0 5 10 15 20 SUM OF SQUARE OF DEVIATIONS FOR THIS ERROR CURVE=3.0490262 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 −10 −5 0 5 10 15 20 SUM OF SQUARE OF DEVIATIONS FOR THIS ERROR CURVE=3.1182013 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 −10 −5 0 5 10 15 20

Fig 12.2

6.4 Solution to ODE where analytical is solution is dicult Example 11: Consider a problem

dy 2 + 18t + 68t2 + 180t3 + 250t4 + 250t5 = dt y2 for initial condition y(t=1)=0, and 0
6.4.1 Solution of Example 11 by Scialb ode solver //**S //##S de('[z]=g(t,y)','z=(2+18*t+68*t^2+180*t^3+250*t^4+250*t^4)/(y^2)') t0=0; y0=1; tn=24; h=0.1; t=0:0.1:24; y=ode(y0,t0,t,g); //We have already dene f(t,y) in previous solution. plot(t,y) xtitle('Solution plot to Example 11','t','y')

51

//##E //**E

Solution plot to Example 11 1400

1200

1000

y

800

600

400

200

0 0

5

10

15

20

25

t

Fig 13

6.5 Modeling and Simulation of a Pressure wave generator (One of the KVPY mentorship holder work on the same project in the Refrigeration Lab at Department of Mechanical Engineering of IIT Bombay during KVPY summer camp 2009. We have model the system and simulate for optimize result).

since there are spring attached (not shown) in the system for restoring purpose and it has found that cylinders also behave like spring, so if we have have combined spring constant as k, then natural frequency of the system

w0 =

p k/m

where m is the mass of the moving part of the system. In this derivation damping forces due to friction and others are assumed to be zero. Since a current carrying coil is suspended in the magnetic eld, so a Lorentz force appear as

52

Fl =

2πrnBi0 (L − x)sin(wt) m

where

r =radius of coil n =no. of turns per unit length of coil B =Magnetic eld produce by Magnetic core i0 =Amplitude of supplied AC current to coil L =Length of the coil inside magnetic core at x =Position of the moving part (or piston) at w =Angular frequency of the AC current

steady or mean state any time

while balancing these force we have the dierential equation as

Bi0 d2 x + w02 x = (L − x)sin(wt) dt2 m or

2πrnBi0 d2 x 2πrnBi0 sin(wt)}x = Lsin(wt) + {w02 + dt2 m m

Let

a = w02 b=

2πrnBi0 m

So we have equations as

d2 x + {a + bsin(wt)}x = bLsin(wt) dt2 Above things are for ideal conditions, in practical we have one more term in the equation as

d2 x b + {a + bsin(wt)}x = (2L + hg )sin(wt) dt2 2 where

hg is

the height of air gap between coil and magnetic core.

To solve this dierential equations in Scilab we can re-write them as

dx1 = x2 dt and

dx2 b = −{a + bsin(wt)}x1 + (2L + hg )sin(wt) dt 2

Following program have been made to solve above dierential equations

53

Program 13: //##S function []=mech1(n,r,B,w,i0,m,v0,t0,tn,Dt,hg,k,w0,L) if length(hu)<>length(w0) then error('length of w0 and lv must be same') abort; end t=t0:Dt:tn; p=length(w0); a=(w0)^2; b=B*(i0)*n*2*%pi*r/m; de('[z]=f(t,x)','z=[x(2);-(a+b*sin(w*t))*x(1)+(b/2)*(hg+2*L)*sin(w*t)]') t0=[0; 0]; x=ode([0 v0]',0,t,f ); xset('window',k) clf subplot(2,1,1) plot(t,x(1,:)) xtitle('Position Vs time plot','time','Position(x)') subplot(2,1,2) plot2d(x(1,:)',x(2,:)') xtitle('Phase plot between velocity and position','Position','Velocity') endfunction ##E In the program v0 is the initial condition for the velocity, t0 and tn are start and end time for solution, Dt is the time increment for time range, k is the gure window no. on which we want the current solution plot. For ideal conditions of manufacturing the solution plot should have sinusoidal pattern, which depends upon all input parameteres. Simulation for the above has been done for following parameters (all values are given in their corresponding SI units) to approach ideal conditions and the corresponding solutions are plotted in Fig 13:

−1

(a)n=12000m

,r=0.02m,B=0.5T,w=314rad/sec,

i0 = 2amp,m=0.03kg,v0 =0.25m/sec,t0 =0sec,tn =25sec, Dt==0.05sec,hg =0.01m,k=0,w0 =250rad/sec,L=0.1m.

−1

(b)n=12000m

,r=0.02m,B=0.5T,w=314rad/sec,

i0 = 2amp,m=0.03kg,v0 =0.25m/sec,t0 =0sec,tn =25sec ,Dt==0.05sec,hg =0.01m,k=1,w0 =350rad/sec,L=0.1m. (c)n=12000m

−1

,r=0.02m,B=0.5T,w=314rad/sec,

i0 = 2amp,m=0.03kg,v0 =0.25m/sec,t0 =0sec,tn =100sec ,Dt==0.05sec,hg =0.01m,k=2,w0 =330rad/sec,L=0.15m.

54

Fig 14 Fig 13 clearly explains how changes in parameters can change the solution.

7 Demerits and diculties in Scilab and suggestion for improvement (Scilab has been compared with MatLab for same problems and tasks. However, it is found that Scilab has some feature better than Matlab but here only demerits and diculties found are shown for its development). 1. It is limited to the developer community and open source users. It needs advertisement and implementation in various disciplines of Engineering and Science for better feedback and development. 2. The implementation of some eld especially studied ODE is quite simple and readymade functions are available. function implementations or when no.

The program crashes after 4-5 of iterations are large.

This has

been especially seen during implementation of a Chemical Engineering ODE of sti type. 3. Scilab has a very poor help browser as compared to its competitive software packages. While searching by some string, mostly Scilab only looks for matching strings in the topic title only and not inside the topic also. The help contains neither a categorical listing of functions nor an extensive search. It does not feature tutorials or a Getting Started tool. The online help available is also limited and not many tutorials are available.

55

4. Scilab has week GUI in command, Editor and other windows while comparing with Matlab. While error for something wrong in Editor window, Scilab only shows at which line it is and not shows the column position, for complicated functions sometimes it becomes dicult to x it. Similar case is with brackets in the Editor window, at any moment it doesn't shows desired animation for paired brackets. 5. Limited no. of toolboxes availability (if available then dicult to access) lacks Scilab. For example a basic toolbox for computation of Laplace and Inverse Laplace Transform which comes frequently in Electrical, Mechanical, Chemical Engineering disciplines. The spectrum of topics dealt with, even in the contributors page is not large. 6. Many times it crash or hang or shows 'stack size exceeded' for large no. of iterations. 7. A good timer function is not available for comparative studies. Example: comparison of Runge-Kutta, Fehlberg's Runge-Kutta, Adams-Bashforth and Sti methods for the solution of ODEs. 8. Scilab don't have symbolic solution environment that means readymade functions like sin, log, sinh etc. which reduces its capacity to solve symbolic equations. Dicult to dene functions with symbolic terms, in the manner we dene polynomial by 'poly' command. Example

−− > y = poly([2, 43, 4], x,0 coef f 0 ) −− > derivat(y) −− > ans = 43 + 8x Now how to derivate this

y = cosx + eax lnx similar thing is for integration. 9. Frequently while loading les it hangs and request for unnecessary time. 10. While looking for help regarding plot or scicos related commands, help browser doesn't show visual graphic example like plots and scicos model.

56

Summary and Conclusion The purpose of the present investigation was to study an optimize approach for solution of a particular category of ordinary dierential equation. However, almost all type of ordinary dierential equation can solve by all given methods but our concern is on two points rst is the simulation time for numerical solution and the second is the time step required for a desired accuracy of solution. In general for a particular method time taken for computation increases as time step (or may be independent variable step in some cases) reduces. Comparative study has done by keeping one of this constant and other as variable.

It is found that although Euler's method can

solve problems but simulation time required by it is usually more or say time step should be less in this case, especially in case of sti category problem it takes large time as compared to other methods. On, other hands in this study methods other than Euler's method all methods are found to more or less similar (Note: Higher Engineering problems are not taken in this study they may distinguish these methods) until unless problem is not of sti type. In the sti category problem it is possible to solve them other methods but simulation time are very high for a desired accuracy which is found to be highest in case of Euler's method, Sti solver stands best here in terms of simulation time, time step and performance of machine. Meanwhile time of above work some other task related to main work done for example curve tting, solution of ode by Laplace Transform Technique which concludes that Scilab generates interactive plots at one place for a better analysis but it should be equipped with some essential toolboxes.

The main conclusion of the approach is that

Scilab is capable to lead everyone but needs a sincere and successive development for this achievement. Also Scilab is capable to solve different kind of ordinary dierential equations and their system which regularly comes in industry.

The important consequence is being

the Scilab free of cost and open source so that can be use easily and sometimes it is faster than other software packages.

57

References 1. http://www.physics.arizona.edu/~restrepo/475B/Notes/source/node16.html, http://en.wikipedia.org/wiki/sti_equation 2. numerical and statistical methods with SCILAB for science and engineering, vol. 1, gilberto e. urroz 3. Youngstown state University, http://www.eng.ysu.edu/~jalam/ engr6924s07/sessions/session27/session27.pdf 4. Stanford University, http://uid.stanford.edu/~nger/teaching/ numerical_methods_02/tutorials/tutorial2.pdf 5. http://planetmath.org/encyclopedia/lorenzEquation.html

58

Kishor Vaigyanik Protsahan Yojana

Jun 11, 2009 - window. //**E: Shows the end of the part of program executed online in Scilab win- dow. //##S: Shows start of the writer's contribution to the program or in the matter taken from ... Microsoft Windows XP Professional version 2002 SP2, Intel(R) Pen- tium(R) D CPU ... 1 Scilab. 8. 1.1 SCILAB HISTORY .

2MB Sizes 1 Downloads 232 Views

Recommend Documents

Yojana Hindi March 2018.pdf
There was a problem previewing this document. Retrying... Download. Connect more apps... Try one of the apps below to open or edit this item. Yojana Hindi ...

Navi Vardhit Pension Yojana-GR.pdf
Navi Vardhit Pension Yojana-GR.pdf. Navi Vardhit Pension Yojana-GR.pdf. Open. Extract. Open with. Sign In. Main menu. Displaying Navi Vardhit Pension ...

Navi Vardhit Pension Yojana-GR.pdf
Navi Vardhit Pension Yojana-GR.pdf. Navi Vardhit Pension Yojana-GR.pdf. Open. Extract. Open with. Sign In. Main menu. Displaying Navi Vardhit Pension ...

Yojana&Kurukshetra-Feb-IASbaba.pdf
village people. India has laid emphasis on national e-governance plan which aims to transform India. into digital empowered society and knowledge economy.

Pradhan Mantri Jeevan Jyoti Bima Yojana - UCO Bank
... BIMA YOJANA – CLAIM FORM. (to be completed by the Claimant & Bank). 1. NAME OF THE SCHEME. : Pradhan Mantri Jeevan Jyoti. Bima Yojana. 2. POLICY NO. : 3. FULL NAME AND ADDRESS. OF THE BANK. : 4. NAME OF THE DECEASED MEMBER : 5. SAVINGS BANK ACC

Pradhan Mantri Vaya Vandana Yojana - 842- Proposal Form.pdf ...
Pradhan Mantri Vaya Vandana Yojana - 842- Proposal Form.pdf. Pradhan Mantri Vaya Vandana Yojana - 842- Proposal Form.pdf. Open. Extract. Open with.

Pradhan Mantri Awas Yojana PM Housing Scheme ...
Ministry of Housing & Urban Poverty Alleviation ... Housing for All (Urban) .... Awas Yojana PM Housing Scheme PMAYPMAY_Housing for all Guidelines.pdf.

Devadasan N, Swarup A. Rashtriya Swasthya Bima Yojana - An ...
2007 Insurance Regulatory and Development Authority. Please reproduce with due permission. Unless explicitly stated, the information and views published in this. Journal may not be construed as those of the Insurance Regulatory. and Development Autho

plan 828: varishtha pension bima yojana -
Aug 13, 2014 - This is an immediate Annuity plan, which provides for immediate pension in consideration of lump sum amount. This is a Government subsidized scheme. The pension payable is guaranteed for the life time with return of purchase price. Gua

vardhit-pension-yojana-gr-1.pdf
Page 1 of 5. TFP!q$qZ__5 YL NFB, SZJFDF\ VFJTL GJL GSSL. SZ[, JlW"T 5[gXG IMHGFP. ( NEW DIFINED ONTRIBUTION. PENSION SCHEME. U]HZFT ;ZSFZ. GF6F\ lJEFU4. ;lRJF,I4 UF\WLGUZP. 9ZFJ S|DF\S o G5GvZ__#vÒVMVF.v!_v5L. TFZLB o !(q#qZ__5. S[gã ;ZSFZGL ;[JFDF

pm-ujjwala-yojana-application-form-hindi.pdf
Sign in. Page. 1. /. 2. Loading… Page 1 of 2. Page 1 of 2. Page 2 of 2. Page 2 of 2. pm-ujjwala-yojana-application-form-hindi.pdf.

PMAY Scheme PDF in Hindi Pradhan Mantri Awas Yojana PM ...
How To Apply For Pradhan Mantri. Awas Yojana PMAY in Hindi : अब तक इसके बारे म पूर जानकार ाa नहं ह, | इस pmay योजना से जुड़ने के िलए.

Study of Rashtriya Swasthya Bima Yojana (RSBY) Health Insurance ...
Study of Rashtriya Swasthya Bima Yojana (RSBY) Health Insurance in India, Gujarat District.pdf. Study of Rashtriya Swasthya Bima Yojana (RSBY) Health ...

pm-ujjwala-yojana-application-form-hindi.pdf
pm-ujjwala-yojana-application-form-hindi.pdf. pm-ujjwala-yojana-application-form-hindi.pdf. Open. Extract. Open with. Sign In. Main menu. Displaying ...

Yojana Decembe 2015 FINAL PDF-c.pdf
Page 3 of 80. YOJANA December 2015 3. Our Representatives : Ahmedabad: Amita Maru, Bengaluru: B.K. Kiranmai, Chennai: A. Elangovan, Guwahati: ...

Evergreen SANSAD ADARSH GRAM YOJANA by Deepak Suchde.pdf
There was a problem previewing this document. Retrying... Download. Connect more apps... Try one of the apps below to open or edit this item. Evergreen ...