Could anyone help me in fixing the problem in my code for solving the initial value problem numerically: y'=y^2, y(0)=1?
2 views (last 30 days)
Show older comments
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
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
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?
Accepted Answer
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)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!