# Could anyone help me in fixing the problem in my code for solving the initial value problem numerically: y'=y^2, y(0)=1?

10 views (last 30 days)
Iqbal Batiha on 29 Jan 2024
Commented: Iqbal Batiha on 29 Jan 2024
clc;
close all;
clear all;
syms w wp wpp wppp
% Define the differential equation: y' = y^2, y(0)=1, t in [0,1)
f = @(t,y) y^2;
fp = @(t,y) 2*y^3;
fpp = @(t,y) 6*y^4;
fppp = @(t,y) 24*y^5;
% Choose the step size h and create the vector
% of t values from t0 to tf incremented by h
t0 = 0;
tf = 0.9;
h = 0.09;
t = t0:h:tf;
% Plot the approximation with the exact solution
t_exact = t0:h:tf;
y_exact = 1./(1-t_exact);
% Initialize a vector of y values as a zero vector
% and set the initial value: y(t0) = y0
y = zeros(size(t));
y0 = 0.5;
y(1) = y0;
for n = 1:(length(t)-1)
k = f(t(n),y(n));
kp = fp(t(n),y(n));
kpp = fpp(t(n),y(n));
kppp = fppp(t(n),y(n));
eqn1 = f( t(n)+h, y(n)+ h*( (1/2)*k + ((3*h)/28)*kp + ((h^2)/84)*kpp + ((h^3)/1680)*kppp ) ...
+ h*( (1/2)*w - ((3*h)/28)*wp + ((h^2)/84)*wpp - ((h^3)/1680)*wppp ) ) == w;
eqn2 = fp( t(n)+h, y(n)+ h*( (1/2)*k + ((3*h)/28)*kp + ((h^2)/84)*kpp + ((h^3)/1680)*kppp ) ...
+ h*( (1/2)*w - ((3*h)/28)*wp + ((h^2)/84)*wpp - ((h^3)/1680)*wppp ) ) == wp;
eqn3 = fpp( t(n)+h, y(n)+ h*( (1/2)*k + ((3*h)/28)*kp + ((h^2)/84)*kpp + ((h^3)/1680)*kppp ) ...
+ h*( (1/2)*w - ((3*h)/28)*wp + ((h^2)/84)*wpp - ((h^3)/1680)*wppp ) ) == wpp;
eqn4 = fppp( t(n)+h, y(n)+ h*( (1/2)*k + ((3*h)/28)*kp + ((h^2)/84)*kpp + ((h^3)/1680)*kppp ) ...
+ h*( (1/2)*w - ((3*h)/28)*wp + ((h^2)/84)*wpp - ((h^3)/1680)*wppp ) ) == wppp;
sol = solve([eqn1, eqn2, eqn3, eqn4], [w, wp, wpp, wppp])
Sol1 = sol.w;
Sol2 = sol.wp;
Sol3 = sol.wpp;
Sol4 = sol.wppp;
y(n+1) = y(n) + h*( (1/2)*k + ((3*h)/28)*kp + ((h^2)/84)*kpp + ((h^3)/1680)*kppp ) ...
+ h*( (1/2)*Sol1 - ((3*h)/28)*Sol2 + ((h^2)/84)*Sol3 - ((h^3)/1680)*Sol4 )
end
sol = struct with fields:
w: [5×1 sym] wp: [5×1 sym] wpp: [5×1 sym] wppp: [5×1 sym]
Unable to perform assignment because the left and right sides have a different number of elements.
y1 = zeros(size(t));
y01 = 0.5;
y1(1) = y01;
for n = 1:(length(t)-1)
k11 = f(t(n),y1(n));
y1(n+1)=y1(n) + h*k11+(h^2/2)*fp(t(n),y1(n))+(h^3/6)*fpp(t(n),y1(n))+(h^4/24)*fppp(t(n),y1(n))
end
figure(1)
plot(t,y,'r*','linewidth',4)
hold on
plot(t,y1,'g','linewidth',4)
hold on
plot(t_exact,y_exact,'b','linewidth',2)
legend('Obreschkoff solution of order 4','Taylor solution of order 4','Exact solution')
grid on
xlabel('t','fontsize',14);
ylabel('y(t)','fontsize',14);
title('Exact vs. Numerical Solutions','fontsize',14);
h1=gca;
set(h1,'fontsize',14);
fh1 = figure(1);
set(fh1, 'color', 'white')
E1 = abs(y_exact - y)
E2 = abs(y_exact - y1)
figure(2)
plot(t,E1,'r','linewidth',2);
hold on
plot(t,E2,'b','linewidth',2);
legend('Obreschkoff method of order 4','Taylor method of order 4')
grid on
xlabel('t','fontsize',14);
ylabel('Error','fontsize',14);
title('Absolute error','fontsize',14);
h2=gca;
set(h2,'fontsize',14);
fh2 = figure(2);
set(fh2, 'color', 'white')
James Tursa on 29 Jan 2024
Th@Iqbal Batiha The title of this post is to solve the problem numerically, but your code makes an attempt at some type of symbolic solution. So, which is it that you want? A numeric or a symbolic solution?
Iqbal Batiha on 29 Jan 2024
Hello Sir! Actually, I'm working on establishing a new method that is better than the Runge-Kutta method. We hope to do well in that thing!

Alan Stevens on 29 Jan 2024
Here's an attempt using fminsearch.
I don't know anything about the Obreschkoff method, so couldn't say if your equations are incorrect or my modifications have messed things up!
Note that your "exact" equation is only true when y(0) = 1, not when y(0) = 0.5.
% Define the differential equation: y' = y^2, y(0)=1, t in [0,1)
f = @(y) y^2;
fp = @(y) 2*y^3;
fpp = @(y) 6*y^4;
fppp = @(y) 24*y^5;
% Choose the step size h and create the vector
% of t values from t0 to tf incremented by h
t0 = 0;
tf = 0.9;
h = 0.09;
t = t0:h:tf;
% Plot the approximation with the exact solution
y_exact = 1./(1-t); % Only true if y(0) = 1
% Initialize a vector of y values as a zero vector
% and set the initial value: y(t0) = y0
y = zeros(size(t));
y0 = 1;
y(1) = y0;
w = [0; 1; 2; 3];
for n = 1:(length(t)-1)
k = f(y(n));
kp = fp(y(n));
kpp = fpp(y(n));
kppp = fppp(y(n));
eqn = @(w) abs(f( y(n)+ h*( (1/2)*k + ((3*h)/28)*kp + ((h^2)/84)*kpp + ((h^3)/1680)*kppp ) ...
+ h*( (1/2)*w(1) - ((3*h)/28)*w(2) + ((h^2)/84)*w(3) - ((h^3)/1680)*w(4) ) ) - w(1)) ...
+abs(fp( y(n)+ h*( (1/2)*k + ((3*h)/28)*kp + ((h^2)/84)*kpp + ((h^3)/1680)*kppp ) ...
+ h*( (1/2)*w(1) - ((3*h)/28)*w(2) + ((h^2)/84)*w(3) - ((h^3)/1680)*w(4) ) ) - w(2)) ...
+abs(fpp( y(n)+ h*( (1/2)*k + ((3*h)/28)*kp + ((h^2)/84)*kpp + ((h^3)/1680)*kppp ) ...
+ h*( (1/2)*w(1) - ((3*h)/28)*w(2) + ((h^2)/84)*w(3) - ((h^3)/1680)*w(4) ) ) - w(3)) ...
+abs(fppp( y(n)+ h*( (1/2)*k + ((3*h)/28)*kp + ((h^2)/84)*kpp + ((h^3)/1680)*kppp ) ...
+ h*( (1/2)*w(1) - ((3*h)/28)*w(2) + ((h^2)/84)*w(3) - ((h^3)/1680)*w(4) ) ) - w(4));
Sol = fminsearch(eqn, w);
y(n+1) = y(n) + h*( (1/2)*k + ((3*h)/28)*kp + ((h^2)/84)*kpp + ((h^3)/1680)*kppp ) ...
+ h*( (1/2)*Sol(1) - ((3*h)/28)*Sol(2) + ((h^2)/84)*Sol(3) - ((h^3)/1680)*Sol(4) );
end
y1 = zeros(size(t));
y01 = 1;
y1(1) = y01;
for n = 1:(length(t)-1)
k11 = f(y1(n));
y1(n+1)=y1(n) + h*k11+(h^2/2)*fp(y1(n))+(h^3/6)*fpp(y1(n))+(h^4/24)*fppp(y1(n));
end
figure(1)
plot(t,y,'r*','linewidth',4)
hold on
plot(t,y1,'bo','linewidth',4)
hold on
plot(t,y_exact,'g','linewidth',2)
legend('Obreschkoff solution of order 4','Taylor solution of order 4','Exact solution')
grid on
xlabel('t','fontsize',14);
ylabel('y(t)','fontsize',14);
title('Exact vs. Numerical Solutions','fontsize',14);
h1=gca;
set(h1,'fontsize',14);
fh1 = figure(1);
set(fh1, 'color', 'white')
E1 = abs(y_exact - y);
E2 = abs(y_exact - y1);
figure(2)
plot(t,E1,'r*','linewidth',2);
hold on
plot(t,E2,'bo','linewidth',2);
legend('Obreschkoff method of order 4','Taylor method of order 4')
grid on
xlabel('t','fontsize',14);
ylabel('Error','fontsize',14);
title('Absolute error','fontsize',14);
h2=gca;
set(h2,'fontsize',14);
fh2 = figure(2);
set(fh2, 'color', 'white')
Iqbal Batiha on 29 Jan 2024
Yap! You are right! I see that the code is working well now! Thanks a lot for your help and consideration!