Download MATLAB Programming for Solving Ordinary Differential Equations using Runge-Kutta Methods and more Study notes Mathematical Methods for Numerical Analysis and Optimization in PDF only on Docsity!
EXPT-1: Introduction to MATLAB programming. Solving physical problems
governed by ODE by RK2 and RK4 using MATLAB programming.
1. Introduction:
MATLAB program is a high level interactive graphics oriented high level software language
widely used in engineering design, modeling, simulation, R&D, etc. This program is having
numerous features like math computation, numerical analysis ,electrical circuit analysis,
Algorithm development, scientific and engineering graphics, data analysis and
visualization, application development, etc. Further, the program is having another
significant interactive built-in ‘SIMULINK’ program which is operated as ‘drag and drop
tool’ for the similar studies.
2. Objective:
In this experiment students will learn the art of MATLAB programming and apply this for
studies and analysis of various numerical optimization techniques Runge-Kutta 2nd^ and 4th
order methods.
3. Requirements of Apparatus:
(i) LAPTOP/Desktop PC
(ii) MATLAB/Simulink Program by MathWorks Inc.
In order to study and learn the MATLAB/SIMULINK program, it is to be downloaded
(whatever latest version is available) in the PC/Labtop of each student. The memory
requirement for the program may be around 15-20 GW. Once loaded and operationalize the
software, students can work on that.
4. Virtual Lab Experiment:
After the program is clicked, there will be prompt to be displayed (>>) on the screen. As the
MATLAB has various built in commands, these can be practiced with the help of ‘help’
command. However, it is to be noted that in MATLAB everything is a matrix, and
MATLAB stands for MATrix LABoratory After being familiar with the commands, the
RK-2 and RK-4 numerical techniques are to be studied by developing MATLAB code
successfully. However, the important features of the MATLAB programming can be learnt
by practicing the followings:
4.1 Typical Command window on the screen:
4.2 File Names
5.1 VARIABLES
➢ Symbolic variable : To create symbolic expressions, first create symbolic variables, and then use operations on them. The syms command is a convenient shorthand for the ‘sym’ syntax. Declare variables x, y, m as symbolic i.e. syms x, y, m. For example:
syms x y E = x^2+6x+7; F = subs(E,x,y) F = y^2+6y+ ➢ MATLAB Pre-defined Constants (built-in) ➢ Statement or Row vector ‘ too long to fit in one line’:
Line can be continued to the next line by typing 3 or more periods ( i.e. ‘. ‘), then pressing < enter>
➢ Commands for conversion from Polar to Cartesian coordinates and vice - versa ➢ Logical operators: AND, OR, NOT ➢ Trigonometric functions : ➢ Differenttiation and integration commands
➢ Polynomial function Commands ➢ Some Built-in Functions :
6. Flow Control Loops/ Loop Commands:
- For-End
- While-End
- If-else
- If-elseif-else-end
- Switch-case-otherwise-end 7. I/O Commands:
- ‘input’ command display text on the screen,waits for the user to enter something from the keyboard and stores in the specified variable e.g X=input(‘Pl enter value of x=’)
- Output commands disp (A), d isp(‘ Hallo, How are You’) [display values or text on the window] fprintf (format), fprintf (format, variables), fprintf (fid, format, variables), [ fprintf displays or prints the output data to the command window or a file]
9.0 M-Files:
When problems become complicated and require reevaluation, entering command at
MATLAB prompt is not practical and therefore,
- In addition to the above, there is another function called anonymous function (discussed
later)
- Do not give script file name as a MATLAB command or function or same variable name it computes
➢ User Defined Function File (.m file) & Function Header
- Syntax of a function function out1= functionname1 (inp1) function out1= functionname2 (inp1,inp2,inp3) function [out1,out2]= functionname 3 (inp1,inp2,inp3,inp4) [This command should be written at the beginning of the m-file and saved with a file name same as the function name (e.g. functionname1 ) ;Output parameter 0 or 1, can omit [ ] ; 2 or more parameters separated by commas(,) ] Examples:
R=input('Enter R: '); C=input('Enter C: '); L=input('Enter L: '); w=input('Enter w: '); y=impedance(R, C, L, w); fprintf ('\n The magnitude of the impedance at %.1f rad/s is %.3f ohm\n', w, y(3)); fprintf ('\n The angle of the impedance at %.1f rad/s is %.3f degrees\n\n', w, y(4)); Coding: % function without brakets- whether the give no. lies between 10 & 100 or not function inrange a=input(‘enter the number=\n’); If ((a>10) & (a<100)) disp(‘The number lies betn 10 & 100’); else disp(‘The number does not lie betn 10 & 100’); end; 9 .1 Anonymous function : Creating simple functions without having to create m-files each time. Anonymous function can be constructed either at the MATLAB command line or in any m-file function or script.
Syntax : fhandle = @ (arg1,arg2) expression
The sign @is the MATLAB operator that constructs a function handle giving a means of invoking the function. Stores the value of the handle in variable fhandle Coding: a=1.3; b=0.2; c=30; parabola=@(x) ax.^2+bx+c; fplot (parabola, [-25 25])**
10. Plot Commands: ( follow any MATLAB reference Book) 11.0 Solving Differential Equations: Syntax: dsolve('eq1','eq2',...,'cond1','cond2',...,'v') eq1, eq2,... are the differential equations wrt the independent variable ,v. Letter ‘D’ denotes differentiation of the function ‘y’ with respect to the independent variable (t). Dy is dy/dt, D2y is d^2 y /dt^2 , D3y is d^3 y /dt^3. ...Dny is dny /dtn^ , indicating first ,second, third and nth derivative of the function (y) wrt the independent variable,t**.
dsolve(’D2y=c^2y’,’y(0)=1’,’Dy(0)=0’)* ans = 1/2exp(ct)+1/2exp(-ct) >>dsolve(’Dx=3x+4y’,’Dy=-4x+3y’, ’x(0)=0’,’y(0)=1’)** [x,y] = x = exp(3t)sin(4t), y = exp(3t)cos(4t) Solving Ordinary Differnential Equation (ODE):
7.0 RK2 and RK4 Methods: The diff equations to be solved in power system stability analysis are nonlinear ordinary differential equations with known initial values: d y /dx = f(x,y) where y is the state vector of n dependent variables and x is the independent variable. Our objective is to solve y as a function of x, with initial values of y and x equal to y 0 and x 0 respectively. There are many algorithms to solve such engineering problems, and Runge-Kutta method is one of them.
7.1 Pre-requisites
Taylor Series: Euler’s Method: ➢ A numerical technique to solve ordinary differential equations of the form
d y /dx = f(x,y), y(0)=0;
Soln.: Solving y at a value of x =xr (=x0+rh) for r=1,2….
Backward Euler’s Method: We start at the initial point (x 0 , y 0 ). Here, the derivative is to be evaluated at point (x + h) instead of at point h. To obtain the next point, (x 1 , y 1 ), we take the derivative at x 1 (not at x 0 !) and extrapolate it at point (x 0 , y 0 ). Hence Backward Euler’s method gives us, y 1 = y 0 + hf(x 1 , y 1 ). Similarly, at (x 1 , y 1 ) we have y 2 = y1 + hf(x 2 , y 2 ). In general, y n+1 = yn + h f(x n+1, yn+1),
The formula giving the value of y is
Where,
k 1 = (slope at the beginning of step size) x h
k 2 = (first approximation to slope at midstep) x h
k 3 = (second approximation to slope at midstep) x h
k 4 = (slop at the end of step) x h
and , incremental value of y = 1/6( k 1 + 2 k 2 +2 k 3 + k 4 ) = weighted average of estimates
based on slopes at the beginning, midpoint and end of the time step.
In general terms,
yn+1 = yn + 1/6( k 1 + 2 k 2 +2 k 3 + k 4 )
where, k 1 = f (xn, yn) * h and so on for k 2 ,k 3 &k 4.
[If f is independent of t, the differential equation is equivalent to a simple integral, then RK4 is Simpson's rule, and If f is independent of t (so called autonomous system, or time-invariant system, especially in physics), and their increments are not computed at all and not passed to function f] 7.3 Example of Runge-Kutta Method Solve
y’ = f ( t, y ) = y − t^2 + 1
y (0) = 0_._ 5 for 0 ≤ t ≤ 2 and h=0. Exact solution for the problem is y = t^2 + 2 t + 1 – ½(et^ )
7.3.1 RK2 :
Step size h = 0****. 5. For t = 0 to t = 1, steps: t 0 = 0, t 1 = 0****. 5, t 2 = 1, t 3 = 1****. 5, t 4 = 2. Step 0: t0 = 0, y0 = 0.
y 1 = y 0 + ( k 1 + k 2 ) / 2 = 1. 375 (RK2)
7.3.2 RK4 :
