i need help on what is wrong in these codes for performing loop

Please am having issue with my codes. The codes supposed to calculate four values for each step size. its only two values out the four that are correct and one value is a constant value throughout and the other is not correct. I do not know what to do. Here is the code and the correct results it suppose display. % Euler approximation for ODE initial value problem % Euler forward method % File prepared by Oladipo Bolaji John. 25/05/2017 % Calculation of h from xinit, xfinal, and n b1 = 0.949; b2 = 3.439; b3 = 18.72; b4 = 37.51; b5 = 1.169;

xmax=9; xinit= 0; yinit= 0.5; %y(xinit)=yinit {y(0)=0} xinit=0; uinit= 0.0; %y(xinit)=uinit {u(0)=0} xinit=0; vinit= 0.0; %v(xinit)=vinit {v(0)=0} xinit=0; winit= 50; %w(xinit)=winit {w(0)=0}

x = xinit:0.5:xmax h=x(2)-x(1); m=length(x);

% for y values y=zeros(1,m); y(1)=yinit; f=@(x,y) (yinit*b1) -((b1*(yinit^2))/b2); for n=1:m-1 y(n+1)= y(n)+h*f(x(n),y(n)) end; disp(y);

% for u values u=zeros(1,m); u(1)=uinit; g=@(x,u) (((yinit*b3*winit)/ (b4+winit))- (0.9082*b5*uinit)); for n=1:m-1 u(n+1)= u(n)+h*g(x(n),u(n)) end; disp(u);

% for v values v=zeros(1,m); v(1)=vinit; k=@(x,v) (b5*uinit); for n=1:m-1 v(n+1)= v(n)+h*k(x(n),v(n)) end; disp(u);

% for w values w=zeros(1,m); w(1)=winit; l=@(x,w) -1.011*((yinit*b3*winit)/( b4+winit)); for n=1:m-1 w(n+1)= w(n)+h*l(x(n),w(n)) end; disp(w);

t y1 y2 y3 y4

     0  0.50000000    0  50
    0.50000000  0.70275596    2.67398011655811  47.2966061021598
    1.00000000  0.90551192    5.34796023311622  44.5932122043195
    1.50000000  1.10826788    8.02194034967432  41.8898183064793
    2.00000000  1.31102384    10.6959204662324  39.1864244086390
    2.50000000  1.51377980    13.3699005827905  36.4830305107988
    3.00000000  1.71653576    16.0438806993486  33.7796366129585
    3.50000000  1.91929172    18.7178608159068  31.0762427151183
    4.00000000  2.12204768    21.3918409324649  28.3728488172780
    4.50000000  2.32480364    24.0658210490230  25.6694549194378
    5.00000000  2.52755961    26.7398011655811  22.9660610215976
    5.50000000  2.73031557    29.4137812821392  20.2626671237573
    6.00000000  2.93307153    32.0877613986973  17.5592732259171
    6.50000000  3.13582749    34.7617415152554  14.8558793280768
    7.00000000  3.33858345    37.4357216318135  12.1524854302366
    7.50000000  3.54133941    40.1097017483716  9.44909153239632
    8.00000000  3.74409537    42.7836818649297  6.74569763455607
    8.50000000  3.94685133    45.4576619814878  4.04230373671583
    9.00000000  4.14960729    48.1316420980459  1.33890983887558

[bolaji oladipo]

this is the second time of asking this question. I need help please

Answer 1

[Are Mjaavatten]

First, your code is hard to read and needs extensive editing if I copy it to the Matlab editor. Add a blank line before the code and indent the code lines for better readability. And start on a new line after comments.

Your inline functions are nominally functions of two variables, but in reality, the second variable is not used. For instance, you use yinit (which is a constant) where I believe you should use the input variable y. I would suggest to write e.g. g as:

 g=@(x,y,u,v,w) y*b3*w/(b4+w)-0.9082*b5*u

Then use a common loop for all four variables. The expression for u(n+1) would then be:

 u(n+1) = u(n) + g(x(n),y(n),u(n),v(n),w(n))

A vector formulation would be even better. See my answer to a question about Euler's explicit method.

[bolaji oladipo]

thanks will get back to you. I am very grateful and really appreciate your response

[bolaji oladipo]

please I really did not get your comment. can u please kindly send the codes with your inputs to me here? I mean how it suppose to look like.

[Are Mjaavatten]

I have been away, and did not see your comment until yesterday. Here is my modified version of your code:

 b1 = 0.949; b2 = 3.439; b3 = 18.72; b4 = 37.51; b5 = 1.169;
 x = (0:0.5:9)';
 h=x(2)-x(1); 
 m=length(x);

 y = zeros(m,1); y(1) = 0.5;
 u = zeros(m,1); u(1) = 0;
 v = zeros(m,1); v(1) = 0;
 w = zeros(m,1); w(1) = 50;

 f = @(x,y,u,v,w) y*b1 -b1*(y^2)/b2;
 g = @(x,y,u,v,w) y*b3*w/(b4+w)- 0.9082*b5*u; 
 k = @(x,y,u,v,w) b5*u;
 l = @(x,y,u,v,w) -1.011*y*b3*w/(b4+w);

 for n=1:m-1 
    y(n+1)= y(n)+h*f(x(n),y(n),u(n),v(n),w(n));
    u(n+1)= u(n)+h*g(x(n),y(n),u(n),v(n),w(n));
    v(n+1)= v(n)+h*k(x(n),y(n),u(n),v(n),w(n));
    w(n+1)= w(n)+h*l(x(n),y(n),u(n),v(n),w(n));
 end; 
 disp([x,y,u,v,w])

Alternatively, and perhaps more elegantly, you may collect y,u,v,w in a matrix, z:

 z = zeros(m,4);
 z0 = [yinit,uinit,vinit,winit];
 z(1,:) = z0;
 f_vector = @(x,z) [b1*(z(1)-(z(1)^2)/b2); ...
    z(1)*b3*z(4)/(b4+z(4))- 0.9082*b5*z(2); ...
    b5*z(2);-1.011*z(1)*b3*z(4)/(b4+z(4))];
 for n=1:m-1 
    z(n+1,:)= z(n,:)+h*f_vector(x(n),z(n,:))';
 end; 
 disp([x,z])

The Euler method gives a rather course solution to the ODE system. Here is how to solve it using the built-in ode45 solver:

 [x_ode,z_ode] = ode45(f_vector,x,z0);
 disp([x_ode,z_ode]);