ASSIGNMENT-1

1.Interpolation(Newton’s Forward ,Newton’s Backward & Lagrange’s )method.

Solution:

%newtons forward differennce method

n=input('number of points= ');
h=input('h= ');
xp=input('the point of interpolation= ');
x=zeros(1,n);
y=zeros(n,n);
for l=1:n
x(l)=input('x= ');
y(l,1)=input('y= ');
end
for j=2:n
for i=1:n-j+1
y(i,j)=y(i+1,j-1)-y(i,j-1)
end
end
sum_f=0;
p=(xp-x(1))/h;
for k=2:n
term=(prod(p-(0:k-2))*y(1,k))/factorial(k-1);
sum_f=sum_f+term;
end
y_xp_f=y(1,1)+sum_f;
fprintf('the value of y(%f) is = %f\n',xp,y_xp_f);

INPUT:

>> Sol_Q1

number of points= 4

h= 2

the point of interpolation= 8

x= 1
y= 24
x= 3
y= 120
x= 5
y= 336
x= 7
y= 720
OutPUT:
y=
24 96 120 48
120 216 168 0
336 384 0 0
720 0 0 0
the value of y(8.000000) is = 990.000000

%newtons backward differennce method

n=input('number of points= ');
h=input('h= ');
xp=input('point of interpolation= ');
x=zeros(1,n);
y=zeros(n,n);
for l=1:n
x(l)=input('x= ');
y(l,1)=input('y= ');
end
for j=2:n
for i=1:n-j+1
y(n+1-i,j)=y(n+1-i,j-1)-y(n-i,j-1);
end
end
p=(xp-x(n))/h;
sum_b=0;
for k=2:n
term=(prod(p+(0:k-2))*y(n,k))/factorial(k-1);
sum_b=sum_b+term;
end
y_xp_b=y(n,1)+sum_b;
disp(y);
fprintf('the value of y(%f) is = %f\n',xp,y_xp_b);

INPUT:

>> sol_Q1_backw
number of points= 4
h= 2
point of interpolation= 8
x= 1
y= 24
x= 3
y= 120
x= 5
y= 336
x= 7
y= 720

OutPUT:

24 0 0 0
120 96 0 0
336 216 120 0
720 384 168 48
the value of y(8.000000) is = 990.000000

%lagrangr's Method
n=input('number of points= ');
xp=input('point of interpolation= ');
for i=1:n
x(i)=input('x= ');
y(i)=input('y= ');
end
sum=0;
for k=1:n
sum=sum+((prod(xp-x(1:n))*y(k))/(prod(x(k)-x(1:k-
1))*prod(x(k)-(x(k+1:n)))*(xp-x(k))));
end
fprintf('the value of y(%f) is = %f\n',xp,sum);
INPUT:
>> Sol_Q1_lagrag
number of points= 4
point of interpolation= 8
x= 1
y= 24
x= 3
y= 120
x= 5
y= 336
x= 7
y= 720
OutPUT:
the value of y(8.000000) is = 990.000000

𝟏.𝟔 𝟏
2.Approximate the integral ∫𝟏.𝟐 (𝒙 + 𝒙) 𝒅𝒙 using
Trapezoidal,Simpson’s1/3,Simpson’s 3/8 & weddle’s rule.

Solution:

%Trapizoidal Method
f=@(x) (x+1/x);
n=input('number of mesh points= ');
x(1)=input('lowe limit= ');
x(n)=input('upper limit= ');
sum=0;
h=(x(n)-x(1))/(n-1);
for i=2:n-1
x(i)=x(i-1)+h;
end
for j=1:n
y(j)=f(x(j));
end
for k=2:n-1
sum=sum+y(k);
end
trapezoidal=(h/2)*(y(1)+(2*sum)+y(n));
fprintf('the value of the given integration=
%f\n',trapezoidal);

INPUT:
number of mesh points= 10
lowe limit= 1.2
upper limit= 1.6
OUTPUT:
the value of the given integration= 0.847732

%simpson's 1/3 rule

f=@(x) (x+1/x);
n=input('number of mesh points= ');
x(1)=input('lowe limit= ');
x(n)=input('upper limit= ');
sum1=0;sum2=0;
h=(x(n)-x(1))/(n-1);
for i=2:n-1
x(i)=x(i-1)+h;
end
for j=1:n
y(j)=f(x(j));
end
for k=2:n-1
if mod(k,2)==0
sum1=sum1+y(k);
else
sum2=sum2+y(k);
end
end
simp13=(h/3)*(y(1)+(4*sum1)+(2*sum2)+y(n));
fprintf('the value of the given integration=
%f\n',simp13);

INPUT:
number of mesh points= 10
lowe limit= 1.2
upper limit= 1.6
OUTPUT:
the value of the given integration= 0.814920
%simpson's 3/8 rule
f=@(x) (x+1/x);
n=input('number of mesh points= ');
x(1)=input('lowe limit= ');
x(n)=input('upper limit= ');
sum1=0;sum2=0;
h=(x(n)-x(1))/(n-1);
for i=2:n-1
x(i)=x(i-1)+h;
end
for j=1:n
y(j)=f(x(j));
end
for k=2:n-1
if mod(k,3)==1
sum1=sum1+y(k);
else
sum2=sum2+y(k);
end
end
simp38=(3*h/8)*(y(1)+(2*sum1)+(3*sum2)+y(n));
fprintf('the value of the given integration=
%f\n',simp38);

INPUT:
number of mesh points= 10
lowe limit= 1.2
upper limit= 1.6
OUTPUT:
the value of the given integration= 0.847682
%weddle's rule
f=@(x) x+1/x;
n=input('number of mesh points= ');
x(1)=input('lowe limit= ');
x(n)=input('upper limit= ');
sum=0;
y=zeros(1000);
h=(x(n)-x(1))/(n-1);
for i=2:n-1
x(i)=x(i-1)+h;
end
for j=1:n-1
y(j)=f(x(j));
end
for k=2:6:n-1

sum=sum+((5*y(k))+y(k+1)+(6*y(k+2))+y(k+3)+(5*y(k+4))+(2*y(k
+5)));
end
y(n)=f(x(n));
wed=((3*h)/10)*(y(1)+sum+y(n))

number of mesh points= 7

lowe limit= 1.2
upper limit= 1.6
OUTPUT:
the value of the given integration= 0.8477

3.Find the Numerical solution of 𝒚′ = 𝒚 − 𝒕𝟐 + 𝟏, 𝟎 ≤ 𝒕 ≤ 𝟐, 𝒚(𝟎) =

𝟎. 𝟓 𝒘𝒊𝒕𝒉 𝒉 = 𝟎. 𝟐 using Taylor’s method Euler’s method, Modified Euler’s
method,R-K method of order 2&4.
Solution:

f=@(t,y) y-t^2+1;
d=@(t,y) y-t^2+1-(2*t); %derivative w.r.to t to find T(2)
for taylor's method
a=input('lower limit of interval= ');
b=input('upper limit of interval= ');
h=input('h= ');
y(1)=input('initial value of y= ');
n=(b-a)/h;
t(1)=a;t(n)=b;
for i=2:n
t(i)=t(i-1)+h;
end
disp('taylor Euler modified eu RK2 Rk4')
%taylor's method
for j=2:n
T(j)=f(t(j-1),y(j-1))+(h/2)*d(t(j-1),y(j-1));
y(j)=y(j-1)+(h*T(j));
tyl(j)=y(j);
end
%Euler's method
for k=2:n
y(k)=y(k-1)+(h*f(t(k-1),y(k-1)));
eu(k)=y(k);
end
%Modified Euler's method
for s=2:n
y(s)=y(s-1)+(h/2)*(f(t(s-1),y(s-1))+f(t(s),y(s-
1)+(h*f(t(s-1),y(s-1)))));
meu(s)=y(s);
end
%Rnge kutta Method of order 2
for p=2:n
p1=f(t(p-1),y(p-1));
p2=f(t(p),y(p-1)+(p1*h));
y(p)=y(p-1)+(h/2)*(p1+p2);
rk2(p)=y(p);
end
%Rnge kutta Method of order 4
for l=2:n
k1=h*f(t(l-1),y(l-1));
k2=h*f(t(l-1)+(h/2),y(l-1)+(k1/2));
k3=h*f(t(l-1)+(h/2),y(l-1)+(k2/2));
k4=h*f(t(l),y(l-1)+k3);
y(l)=y(l-1)+(1/6)*(k1+(2*k2)+(2*k3)+k4);
rk4(l)=y(l);
end
for c=2:n
fprintf( ' %f %f %f %f %f
\n',tyl(c),eu(c),meu(c),rk2(c),rk4(c));
end
INPUT:
lower limit of interval= 0
upper limit of interval= 2
h= 0.5
initial value of y= 0.5
OUTPUT:
taylor Euler modified eu RK2 Rk4
1.437500 1.250000 1.375000 1.375000 1.425130
2.679688 2.250000 2.515625 2.515625 2.639603
4.104492 3.375000 3.775391 3.775391 4.006819

4. Solve the following initial value problem over the interval [0,2].Display all your
results graphically.
(a) ode solver .
(b)Using Euler’s method with h=0.5 & h= 0.25
(c) Using Midpoint method with h=0.5 .
(d) Using the fourth order R-K method with h=0.5 .
Solution:
f=@(t,y) (y*t^3)-(1.5*y);
a=input('lower limit of interval= ');
b=input('upper limit of interval= ');
h=input('h= ');
y(1)=input('initial value of y= ');
n=(b-a)/h;
t(1)=a;t(n+1)=b;
for r=2:n
t(r)=a+(r*h);
end
eu(1)=y(1);
mp(1)=y(1);
rk4(1)=y(1);
%(a) using ODE Solver
ode45(f,t,y(1))
%(b) Euler's method
hold on
for i=2:n+1
y(i)=y(i-1)+(h*f(t(i-1),y(i-1)));
eu(i)=y(i);
end
%(c) Mid point method
for j=2:n+1
y(j)=y(j-1)+h*f(t(j-1)+(h/2),y(j-1)+(h/2)*f(t(j-1),y(j-1)));
mp(j)=y(j);
end
for l=2:n+1
k1=h*f(t(l-1),y(l-1));
k2=h*f(t(l-1)+(h/2),y(l-1)+(k1/2));
k3=h*f(t(l-1)+(h/2),y(l-1)+(k2/2));
k4=h*f(t(l),y(l-1)+k3);
y(l)=y(l-1)+(1/6)*(k1+(2*k2)+(2*k3)+k4);
rk4(l)=y(l);
end
for c=1:n+1
fprintf('%f %f %f %f \n',t(c),
eu(c),mp(c),rk4(c));
end
plot(t,eu,'k-',t,mp,'b--',t,rk4,'r*:');
INPUT:
lower limit of interval= 0
upper limit of interval= 2
h= 0.25
initial value of y= 1
OUTPUT:
0.000000 1.000000 1.000000 1.000000
0.500000 0.625000 0.695709 0.691126
0.750000 0.410156 0.514823 0.506139
1.000000 0.299606 0.422385 0.412689
1.250000 0.262156 0.414844 0.406683
1.500000 0.291853 0.535346 0.538192
1.750000 0.428659 0.996433 1.086770
2.000000 0.842248 2.876751 3.863870
2.000000 2.210900 13.429725 24.072147

𝒅𝟐 𝒅𝒚
5.Solve the differendial equation 𝟐 + 𝟑 + 𝟐𝒚 = 𝟎,subject to the initial
𝒅𝒙 𝒅𝒙
conditions that 𝒚(𝟎) = 𝟏 &𝒚′(𝟎) = 𝟎, from 0 to 1.
Solution:
f=@(x,Y) [Y(2);-3*Y(2)-2*Y(1)];
y0=[1;0];
xspan=[0 1];
[x,Y]=ode45(f,xspan,y0);
disp([x,Y(:,1)])
plot(x,Y(:,1))
INPUT-OUTPUT:
0 1.0000
0.1000 0.9909
0.2000 0.9671
0.3000 0.9328
0.4000 0.8913
0.5000 0.8452
0.6000 0.7964
0.7000 0.7466
0.8000 0.6968
0.9000 0.6478
1.0000 0.6004

6.The population of X from the year 1930 to thr year 2020 is given in the
following table :
Year 1930 1940 1950 1960 1970 1980 1990 2000 2010

Population 249 277 316 350 431 539 689 833 1014
In million

(a) Fit the data with a second -order Ploynomial .Make a plot of the points &
the polynomial.
(b) Fit the data with linear & spline interpolations.Estimate the population of
1995 with linear & spline interpolations.

Solution(a):

year=[1930 1940 1950 1960 1970 1980 1990 2000 2010];

population=[249 277 316 350 431 539 689 833 1014];
ploynomial=polyfit(year,population,2);
value=polyval(ploynomial,year);
plot(year,population,'r*-',year,value,'b^:');
legend('points','ploynomial');
OUT-PUT:

Solution(b):

yr=[1930 1940 1950 1960 1970 1980 1990 2000 2010];

pop=[249 277 316 350 431 539 689 833 1014];
interp_yr=input('enter the year = ');
year=1930:10:2010;
fit_lin=interp1(yr,pop,year,"linear");
fit_spl=interp1(yr,pop,year,'spline');
pop_lin=interp1(yr,pop,interp_yr,'linear');
pop_spl=interp1(yr,pop ,interp_yr,'spline');
fprintf('the poputalin in the year %d is= %f
millions\n',interp_yr,pop_lin);
fprintf('the poputalin in the year %d is= %f
millions\n',interp_yr,pop_spl);
INPUT-OUTPUT:
enter the year = 1995
the poputalin in the year 1995 is= 761.000000 millions
the poputalin in the year 1995 is= 760.939744 millions
 ASSIGNMENT-2
1.Consider the following linear system of equations:

𝟏𝟎𝒙𝟏 − 𝒙𝟐 + 𝟐𝒙𝟑 =𝟔
−𝒙𝟏 + 𝟏𝟏𝒙𝟐 − 𝒙𝟑 + 𝟑𝒙𝟒 = 𝟐𝟓
𝟐𝒙𝟏 − 𝒙𝟐 + 𝟏𝟎𝒙𝟑 − 𝒙𝟒 = −𝟏𝟏
𝟑𝒙𝟐 − 𝒙𝟑 + 𝟖𝒙𝟒 = 𝟏𝟓

Solve the above system ,correct up to 5 decimal places ,with initial guess
𝒙𝟎 = (𝟎, 𝟎, 𝟎, 𝟎) using:
(a) Jacobi Iterative method.
(b) Gauss-seidel iterative method.
(c) SOR iterative method with 𝝎 = 𝟏. 𝟏 .
Solution(a):

%jacobi Iterative method

A=input('enter the coefficient matrix A= ');
b=input('enter the matrix by value ox b= ');
max_it=input('enter the maximum iteration= ');
n=input('enter the order of the matrix n = ');
c=0;
x0=input('enter the initial solution x0=');
x1=zeros(1,n);
tol=input('enter the tolarance= ');
k=0;
disp(' x1 x2 x3 x4');
while k<max_it
k=k+1;
for i=1:n
for j=1:n
if i~=j
c=c+(A(i,j)*x0(j));
end
end
x1(i)=(b(i)-c)/(A(i,i));
c=c-c;
end
disp(x1);
err=abs(x1-x0);
if err<tol
disp('the required root')
disp(x1');
break
else
x0=x1;
end
end

INPUT:

Question_1_a
enter the coefficient matrix A= [10 -1 2 0;-1 11 -1 3;2 -1 10 -1;0 3 -1 8]
enter the matrix by value ox b= [6 25 -11 15]
enter the maximum iteration= 30
enter the order of the matrix n = 4
enter the initial solution x0=[0 0 0 0]
enter the tolarance= 0.00001
OUTPUT:
x1 x2 x3 x4
0.6000 2.2727 -1.1000 1.8750

1.0473 1.7159 -0.8052 0.8852

0.9326 2.0533 -1.0493 1.1309

1.0152 1.9537 -0.9681 0.9738

0.9890 2.0114 -1.0103 1.0214

1.0032 1.9922 -0.9945 0.9944

0.9981 2.0023 -1.0020 1.0036

 1.0006 1.9987 -0.9990 0.9989

0.9997 2.0004 -1.0004 1.0006

1.0001 1.9998 -0.9998 0.9998

0.9999 2.0001 -1.0001 1.0001

1.0000 2.0000 -1.0000 1.0000

1.0000 2.0000 -1.0000 1.0000

1.0000 2.0000 -1.0000 1.0000

1.0000 2.0000 -1.0000 1.0000

1.0000 2.0000 -1.0000 1.0000

the required root

1.0000
2.0000
-1.0000
1.0000
Solution(b):

%Gauss-Seidel Iterative method

A=input('enter the coefficient matrix A= ');
b=input('enter the matrix by value ox b= ');
max_it=input('enter the maximum iteration= ');
n=input('enter the order of the matrix = ');
c=0;
x0=input('enter the initial solution x0=');
x1=zeros(1,n);
tol=input('enter the tolarance= ');
k=0;
disp('x1 x2 x3 x4');
while k<max_it
k=k+1;
for i=1:n
for j=1:n
if i~=j
if i<j
c=c+(A(i,j)*x0(j));
end
if i>j
c=c+(A(i,j)*x1(j));
end
end
end
x1(i)=(b(i)-c)/(A(i,i));
c=c-c;
end
disp(x1);
err=abs(x1-x0);
if err<tol
disp('the required root')
disp(x1');
break
else
x0=x1;
end
end

INPUT:

enter the coefficient matrix A= [10 -1 2 0;-1 11 -1 3;2 -1 10 -1;0 3 -1 8]

enter the matrix by value ox b= [6 25 -11 15]
enter the maximum iteration= 30
enter the order of the matrix n = 4
enter the initial solution x0=[0 0 0 0]
enter the tolarance= 0.00001
OUTPUT:

x1 x2 x3 x4
0.6000 2.3273 -0.9873 0.8789

1.0302 2.0369 -1.0145 0.9843

1.0066 2.0036 -1.0025 0.9984

1.0009 2.0003 -1.0003 0.9998

1.0001 2.0000 -1.0000 1.0000

1.0000 2.0000 -1.0000 1.0000

1.0000 2.0000 -1.0000 1.0000

the required root

1.0000
2.0000
-1.0000
1.0000
Solution(c):
%SOR method
A=input('enter the coefficient matrix A= ');
b=input('enter the matrix by value ox b= ');
max_it=input('enter the maximum iteration= ');
n=size(A,1);
x0=input('enter the initial solution x0=');
tol=input('enter the tolarance= ');
w=input('enter the relaxation factor w= ');
x1=zeros(1,n);
k=0;
c=0;
disp('x1 x2 x3 x4');
while k<max_it
k=k+1;
for i=1:n
for j=1:n
if i~=j
if i<j
c=c+(A(i,j)*x0(j));
end
if i>j
c=c+(A(i,j)*x1(j));
end
end
end
x1(i)=((1-w)*x0(i))+((w*(b(i)-c))/(A(i,i)));
c=c-c;
end
disp(x1);
err=abs(x1-x0);
if err<tol
disp('the required root')
disp(x1');
break
else
x0=x1;
end
end
fprintf('total number of iteration=%d\n',k);
INPUT
enter the coefficient matrix A= [10 -1 2 0;-1 11 -1 3;2 -1 10 -1;0 3 -1 8]
enter the matrix by value ox b= [6 25 -11 15]
enter the maximum iteration= 30
enter the initial solution x0=[0 0 0 0]
enter the tolarance= 0.00001
enter the relaxation factor w= 1.1
Output:
x1 x2 x3 x4
0.6600 2.5660 -1.0729 0.8565

1.1123 1.9904 -1.0343 1.0136

 0.9952 1.9930 -0.9948 1.0023

0.9986 2.0004 -0.9999 0.9996

1.0002 2.0001 -1.0001 1.0000

1.0000 2.0000 -1.0000 1.0000

1.0000 2.0000 -1.0000 1.0000

1.0000 2.0000 -1.0000 1.0000

the required root
1.0000
2.0000
-1.0000
1.0000
total number of iteration=8

2.consider the system of equations given below:

𝟐𝒙𝟏 − 𝟑𝒙𝟐 + 𝟐𝒙𝟑 = 𝟓
−𝟒𝒙𝟏 + 𝟐𝒙𝟐 − 𝟔𝒙𝟑 = 𝟏𝟒
𝟐𝒙𝟏 + 𝟐𝒙𝟐 + 𝟒𝒙𝟑 = 𝟖
(a) Solve the system using Gaussian Ellimination method.
(b) Solve the system using Gauss-Jordan Ellimination method.

Solution(a):

A=input('enter the Augmented matrix A= ');

n=size(A,1);
% Making diagonal element non-zero
for i=1:n-1
for j=2:n
if A(i,i)==0
t=A(i,:);
A(i,:)=A(j,:);
A(j,:)=t;
end
 end
%Making of upper triamghle matrix(row operation)
for k=i+1:n
A(k,:)=A(k,:)-((A(i,:)*A(k,i))/A(i,i));
end
end
disp(A)
x=zeros(1,n);
c=0;
for l=n:-1:1
for s=1:n
if s~=l
c=c+A(l,s)*x(s);
end
end
x(l)=(1/A(l,l))*(A(l,n+1)-c);
c=c-c;
end
disp(x');
INPUT-OUTPUT:
enter the Augmented matrix A= [2 -3 2 5;-4 2 -6 14;2 2 4 8]

2.0000 -3.0000 2.0000 5.0000

0 -4.0000 -2.0000 24.0000
0 0 -0.5000 33.0000

X1= 109
X2= 27
X3=-66
Solution(b):

A=input('enter the Augmented matrix A= ');

n=size(A,1);
x=zeros(1,n);
% Making diagonal element non-zero
for i=1:n-1
for j=2:n
if A(i,i)==0
t=A(i,:);
 A(i,:)=A(j,:);
A(j,:)=t;
end
end
%Making of upper triamghle matrix(row opperation)
for k=i+1:n
A(k,:)=A(k,:)-((A(i,:)*A(k,i))/A(i,i));
end
end
%making of reduced row-echelon form
for a=n:-1:1
for b=a-1:-1:1
A(b,:)=A(b,:)-(A(a,:)*(A(b,a))/A(a,a));
end
end
for p=1:n
A(p,:)=A(p,:)/A(p,p);
x(p)=A(p,n+1);
end
disp(A)
disp(x')
INPUT-OUTPUT:
enter the Augmented matrix A= [2 -3 2 5;-4 2 -6 14;2 2 4 8]

1 0 0 109
0 1 0 27
0 0 1 -66
X1= 109
X2= 27
X3=-66
3.Use power method to approximate the dominant eigenvalue of the matrix
−𝟒 𝟏𝟒 𝟎
𝑨 = (−𝟓 𝟏𝟑 𝟎 ) ,Let𝒙𝟎 = (𝟏, 𝟏, 𝟏)𝑻
−𝟏 𝟎 𝟐
Solution:

%power method to find Dominant Eigenvalue

A=input('enter the Given matrix A= ');
x0=input('enter the initial vectoer in column:');
n=size(A);
max_it=input('Maximum iteration=');
tol=input('tolarance= ');
l=zeros(1,max_it+1);
for i=1:max_it
D=A*x0;
m=max(D(:));
x1=D/m;
tr=x1';
l(i+1)=(tr*(A*x1))/(tr*x1);
err=abs(l(i)-l(i+1));
if err<tol
fprintf('required eigenvalue is= %f\n',l(i+1));
fprintf('total iterations=%d\n',i);
break
else
x0=x1;
end
end

INPUT:
enter the Given matrix A= [-4 14 0;-5 13 0;-1 0 2]
enter the initial vectoer in column: [1;1;1]
Maximum iteration=30
tolarance= 0.0001
OUTPUT:
required eigenvalue is= 6.000053
total iterations=15

4. Solve the system if ODEs

𝒅𝒙
= 𝒂𝒙 = 𝒃𝒙𝒚
𝒅𝒕
𝒅𝒚
= −𝒄𝒚 + 𝒅𝒙𝒚
𝒅𝒕

Taking 𝒂 = 𝟏. 𝟐, 𝒃 = 𝟎. 𝟔, 𝒄 = 𝟎. 𝟖 & 𝒅 = 𝟎. 𝟑 .

Solution:
a=1.2;
b=0.6;
c=0.8;
d=0.3;
f=@(t,Y) [a*Y(1)-b*Y(1)*Y(2);-c*Y(2)+d*Y(1)*Y(2)];
%Y(1)=x;Y(2)=y
tspan=0:0.1:2;
x0=[2;1];
[t,Y]=ode45(f,tspan,x0);
disp([t,Y(:,1)])

OUTPUT:
0 2.0000
0.1000 2.1249
0.2000 2.2597
0.3000 2.4048
0.4000 2.5603
0.5000 2.7264
0.6000 2.9030
0.7000 3.0897
0.8000 3.2861
0.9000 3.4910
1.0000 3.7033
1.1000 3.9208
1.2000 4.1412
1.3000 4.3610
1.4000 4.5762
1.5000 4.7816
1.6000 4.9715
1.7000 5.1390
1.8000 5.2768
1.9000 5.3771
2.0000 5.4327
5.The van der pl equation is a second order ODE
𝒅𝟐 𝒚 𝟐)
𝒅𝒚
− 𝝁(𝟏 − 𝒚 + 𝒚 = 𝟎, 𝒚(𝟎) = 𝟐, 𝒚′ (𝟎) = 𝟎
𝒅𝒙𝟐 𝒅𝒙
Where 𝝁 > 𝟎 is a scalar parameter .Solve the above euation using shooting
method (ode45 solver) and the shoe the result graphically.
Solution:
mu=1.2;
f=@(x,Y) [Y(2);mu*(1-Y(1)^2)*Y(2)-Y(1)];
xspan=0:0.1:20;
y0=[2;0];
[x,Y]=ode45(f,xspan,y0);
disp([x,Y(:,1)])
plot(x,Y(:,1))
xlabel('x-axis');
ylabel('y-axis');
grid on
OUTPUT:
0 2.0000
0.5000 1.8491
1.0000 1.5595
1.5000 1.1618
2.0000 0.5581
2.5000 -0.5095
3.0000 -1.7657
3.5000 -2.0112
4.0000 -1.8179
4.5000 -1.5125
5.0000 -1.0941
5.5000 -0.4449
6.0000 0.7064
6.5000 1.8571
7.0000 1.9860
7.5000 1.7692
8.0000 1.4480
8.5000 1.0014
9.0000 0.2862
9.5000 -0.9616
10.0000 -1.9397
10.5000 -1.9583
11.0000 -1.7181
11.5000 -1.3789
12.0000 -0.8992
12.5000 -0.1059
13.0000 1.2157
13.5000 1.9880
14.0000 1.9243
14.5000 1.6642
15.0000 1.3057
15.5000 0.7875
16.0000 -0.0954
16.5000 -1.4469
17.0000 -2.0094
17.5000 -1.8854
18.0000 -1.6075
18.5000 -1.2279
19.0000 -0.6650
19.5000 0.3182
20.0000 1.6416