Obtain the numerical solution of the first order ODE dy/dx=y^2*(1-2*x) where y(1)=1 with h=0.1.
find the
solution at points x=0.1,0.2,0.3 using R-K method and find solutions at all other points up to x=3 using milne's
predictor-corrector formula and plot to compare with the analytical solution.
x0=0;
xn=3;
y0=1;
h=0.1;
f=@(x,y) y^2*(1-2*x);
x=x0:h:xn;
n=length(x)-1;
y=zeros(1,n+1);
y(1)=y0;
for i=1:3
k1=h*f(x(i),y(i));
k2=h*f(x(i)+h/2,y(i)+k1/2);
k3=h*f(x(i)+h/2,y(i)
+k2/2);
k4=h*f(x(i)+h,y(i)+k3);
y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4);
end
for i=4:n
y(i+1)=y(i-3)+(4*h/3)*(2*f(x(i),y(i))-(f(x(i-1),y(i-1)))
+2*(f(x(i-2),y(i-2))));
y(i+1)=y(i-1)+(h/3)*(f(x(i+1),y(i+1))+4*f(x(i),y(i))+f(x(i-1),y(i-1)));
end
y4=y(5)
y4 =
1.3158
disp(y);
1.0000 1.0989 1.1905 1.2658 1.3158 1.3334 1.3158 1.2658 1.1904 1.0988 0.
ye=eval(dsolve('Dy=y^2*(1-2*x)','y(0)=1','x'));
Warning: Support for character vector or string inputs will be removed in a future release.
Instead, use syms to declare variables and replace inputs such as dsolve('Dy = -3*y') with syms
y(t); dsolve(diff(y,t) == -3*y).
plot(x,y,'*',x,ye,'r')
grid on
xlabel('x')
ylabel('y')
legend('numerical solution','analytical solution')
1
x0=0;
xn=3;
y0=0;
h=0.2;
f=@(x,y) x-y^2;
x=x0:h:xn;
n=length(x)-1;
y=zeros(1,n+1);
y(1)=y0;
for i=1:3
k1=h*f(x(i),y(i));
k2=h*f(x(i)+h/2,y(i)+k1/2);
k3=h*f(x(i)+h/2,y(i)
+k2/2);
k4=h*f(x(i)+h,y(i)+k3);
y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4);
end
for i=4:n
y(i+1)=y(i-3)+(4*h/3)*(2*f(x(i),y(i))-(f(x(i-1),y(i-1)))
+2*(f(x(i-2),y(i-2))));
y(i+1)=y(i-1)+(h/3)*(f(x(i+1),y(i+1))+4*f(x(i),y(i))+f(x(i-1),y(i-1)));
end
y4=y(5)
2
y4 =
0.3046
disp(y);
0 0.0200 0.0795 0.1762 0.3046 0.4556 0.6177 0.7796 0.9323 1.0709 1.
ye=eval(dsolve('Dy=x-y^2','y(0)=0','x'));
Warning: Support for character vector or string inputs will be removed in a future release.
Instead, use syms to declare variables and replace inputs such as dsolve('Dy = -3*y') with syms
y(t); dsolve(diff(y,t) == -3*y).
plot(x,y,'*',x,ye,'r')
Warning: Imaginary parts of complex X and/or Y arguments ignored.
grid on
xlabel('x')
ylabel('y')
legend('numerical solution','analytical solution')
x0=0;
xn=3;
y0=2;
h=0.1;
f=@(x,y) x^2+(y/2);
3
x=x0:h:xn;
n=length(x)-1;
y=zeros(1,n+1);
y(1)=y0;
for i=1:3
k1=h*f(x(i),y(i));
k2=h*f(x(i)+h/2,y(i)+k1/2);
k3=h*f(x(i)+h/2,y(i)
+k2/2);
k4=h*f(x(i)+h,y(i)+k3);
y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4);
end
for i=4:n
y(i+1)=y(i-3)+(4*h/3)*(2*f(x(i),y(i))-(f(x(i-1),y(i-1)))
+2*(f(x(i-2),y(i-2))));
y(i+1)=y(i-1)+(h/3)*(f(x(i+1),y(i+1))+4*f(x(i),y(i))+f(x(i-1),y(i-1)));
end
y4=y(5)
y4 =
2.4652
disp(y);
2.0000 2.1029 2.2131 2.3330 2.4652 2.6125 2.7775 2.9632 3.1728 3.4096 3.
ye=eval(dsolve('Dy=x^2+(y/2)','y(0)=2','x'));
Warning: Support for character vector or string inputs will be removed in a future release.
Instead, use syms to declare variables and replace inputs such as dsolve('Dy = -3*y') with syms
y(t); dsolve(diff(y,t) == -3*y).
plot(x,y,'*',x,ye,'r')
grid on
xlabel('x')
ylabel('y')
legend('numerical solution','analytical solution')
4
x0=0;
xn=3;
y0=2;
h=0.1;
f=@(x,y) 2*exp(x)-y;
x=x0:h:xn;
n=length(x)-1;
y=zeros(1,n+1);
y(1)=y0;
for i=1:3
k1=h*f(x(i),y(i));
k2=h*f(x(i)+h/2,y(i)+k1/2);
k3=h*f(x(i)+h/2,y(i)
+k2/2);
k4=h*f(x(i)+h,y(i)+k3);
y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4);
end
for i=4:n
y(i+1)=y(i-3)+(4*h/3)*(2*f(x(i),y(i))-(f(x(i-1),y(i-1)))
+2*(f(x(i-2),y(i-2))));
y(i+1)=y(i-1)+(h/3)*(f(x(i+1),y(i+1))+4*f(x(i),y(i))+f(x(i-1),y(i-1)));
end
y4=y(5)
5
y4 =
2.1621
disp(y);
2.0000 2.0100 2.0401 2.0907 2.1621 2.2553 2.3709 2.5103 2.6749 2.8662 3.
ye=eval(dsolve('Dy=2*exp(x)-y','y(0)=2','x'));
Warning: Support for character vector or string inputs will be removed in a future release.
Instead, use syms to declare variables and replace inputs such as dsolve('Dy = -3*y') with syms
y(t); dsolve(diff(y,t) == -3*y).
plot(x,y,'*',x,ye,'r')
grid on
xlabel('x')
ylabel('y')
legend('numerical solution','analytical solution')
