midpoint method to solve an ode

I am using a midpoint method to solve an ode !! I am not sure if my code is right but it gives me 0 errors than it shows me only one point in the graph and Inf on the second column !! i cant understand why is that!!

f = @(t,y)(y/t+1 + 5*((t+1)/(1+25*t^2)));

a = input('Enter left end ponit, a: ');

b = input('Enter right end point, b: ');

n = input('Enter no. of subintervals, n: ');

alpha = input('Enter the initial condition, alpha: ');

h = (b-a)/n;

t=[a zeros(1,n)];

w=[alpha zeros(1,n)];

for i = 1:n+1

t(i+1)=t(i)+h;

wprime=w(i)+(h/2)*f(t(i),w(i));

%w(i+1)=w(i)+(h/2)*f( t(i), w(i)+ f(t(i+1),wprime) );

w(i+1)=w(i)+h*f( t(i)+(h/2), wprime );

fprintf('%5.4f %11.8f\n', t(i), w(i));

plot(t(i),w(i),'r*'); grid on;

xlabel('t values'); ylabel('w values');

hold on;

end

0.0000 1.00000000

0.1000 Inf

0.2000 Inf

0.3000 Inf

0.4000 Inf

0.5000 Inf

0.6000 Inf

0.7000 Inf

0.8000 Inf

0.9000 Inf

1.0000 Inf

[Torsten]

Didn't you learn from your other thread that the brackets in your function definition are set incorrectly ?

Answer 1

[ag]

Hi Anas,

I understand that you are getting “Inf” values in second column output of “f” function, and while trying to plot you are observing just one point in the graph.

• To avoid the “Inf” value in the second column, you can try updating the formula for the function f. For example, below formula doesn’t produce “Inf” values and works fine for me,

  f = @(t,y) (y/(t+1))+ 5*(t+1)/(1+25*t^2); 

• As the plot function is used inside the for loop, it creates a new figure handle each time the loop runs. To avoid that, you should declare a figure handle before the “for loop” starts.

The full working code is shown below:

f = @(t,y)(y/(t+1) + 5*((t+1)/(1+25*t^2)));  
a = 10; %input('Enter left end ponit, a: '); 
b = 20; %input('Enter right end point, b: '); 
n = 100; %input('Enter no. of subintervals, n: ');  
alpha = 0.01; %input('Enter the initial condition, alpha: '); 
h = (b-a)/n;  
t=[a zeros(1,n)]; 
w=[alpha zeros(1,n)]; 
figure; 
hold on; 
for i = 1:n+1 
    t(i+1)=t(i)+h; 
    wprime=w(i)+(h/2)*f(t(i),w(i)); 
    %w(i+1)=w(i)+(h/2)*f( t(i), w(i)+ f(t(i+1),wprime) ); 
    w(i+1)=w(i)+h*f( t(i)+(h/2), wprime ); 
    %fprintf('%5.4f %11.8f\n', t(i), w(i)); 
    plot(t(i),w(i),'r*'); grid on; 
    xlabel('t values'); ylabel('w values'); 
end; 
hold off; 

I’ve assumed values of a, b, n and alpha as 10, 20, 100 and 0.01 respectively, to demonstrate the output of the code.

For more details, please refer to the following documentation: https://www.mathworks.com/help/matlab/ref/plot.html

Hope this helps!

Best Regards,

Aryan Gupta