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

Question

Iqbal Batiha 2024-1-29

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2075516-could-anyone-help-me-in-fixing-the-problem-in-my-code-for-solving-the-initial-value-problem-numerica

评论： Iqbal Batiha 2024-1-29

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 个评论
显示 2更早的评论隐藏 2更早的评论

James Tursa 2024-1-29

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 2024-1-29

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 2024-1-29

在 MATLAB Online 中打开

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Iqbal Batiha 2024-1-29

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

请先登录，再进行评论。

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

WeChat

4 个评论
显示 2更早的评论隐藏 2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论