Clear Filters
Clear Filters

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)
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')
  4 Comments
James Tursa
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
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!

Sign in to comment.

Accepted Answer

Alan Stevens
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')

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!