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

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

[Dyuman Joshi]

There are 5 solution pairs for the equation you have solved (see above) i.e. each of the solutions are a 5 element vector (a 5x1 symbolic vector to be precise).

So, in the command where you update y, there is mismatch rhs is 1x1 whereas lhs is 5x1, thus leading to the error you get.

What can you do?

Well, you could choose any one of the solution to proceed with.

There is not enough information for me to comment more on this.

[Sam Chak]

@Iqbal Batiha, I was wondering what underlying numerical method is being used to solve your initial value problem? Is it some kind of advanced Runge-Kutta method?

I'm curious because the choice of numerical method can have a significant impact on the accuracy and efficiency of the solution.

[James Tursa]

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]

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!

Answer 1

[Alan Stevens]

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]

Yap! You are right! I see that the code is working well now! Thanks a lot for your help and consideration!