midpoint rule for matlab

Question

Hello I was asked to create a matlab code for the midpoint rule. What I have is the code for eulers method, so I have to make some modifications, but I am struggling to do it I have the following

function H = heun(f,a,b,ya,M)
h = (b-a)/M;
T = zeros(1,M+1);
Y = zeros(1,M+1);
T = a:h:b;
Y(1) = ya;
for j = 1 : M
k1=feval(f,T(j),Y(j));
k2=feval(f,T(j+1),Y(j)+h*k1);
Y(j+1)=Y(j)+(h/2)*(k1+k2);
end
H = [T' Y'];
function f = dif1(t,y)
f=(t-y)/2;

So yea my function is y'=(t-y)/2 and the interval and other stuff i define in the command window..

If it is possible to make the for loop into a mid point rule I think that is the way to go, any help would be appreciated.

Is it possible? Yes. However I assume you have other questions, but you need to phrase them for us to know what they are! — Ander Biguri, Oct 04 '18 at 09:07
This question is currently "please do my homework, here is some related code". Please try to implement the correct code yourself and tell us where you're stuck what mathematics you're trying to use, what the specific issue is, including a [mcve] and expected outputs. — Wolfie, Oct 04 '18 at 09:15
I am new to matlab, and I wasn't asking you to solve my entire problem I was asking for help, because I spent 5 hours yesterday to understand and implement the above code.. — Chengdu, Oct 04 '18 at 09:31
Okay, so I now tried to compare my code. for the euler method I translate the code from wi+1 = wi+hf(ti,wi) , and that becomes Y(j+1) = Y(j) + h * feval (f,t(j),Y(j));................. Now for midpoint it shows that yn+1 = yn + hf ( tn + h/2, yn + h/2 f(tn,yn)) and can I translate it just as I did with the eulers method? — Chengdu, Oct 04 '18 at 11:30
What you have you correctly identified as Heun's method, or the explicit trapezoidal rule, or the modified Euler method. The explicit midpoint method is sometimes also called RK2 or improved Euler method. Please clarify what exactly you want. Without modifier I would assume that "midpoint method" means the implicit variant, which is still a slightly different method. — Lutz Lehmann, Oct 04 '18 at 14:19
Thank you for your kind reply. In terms of mathematics the "mid point method/rule" or "modified Euler" is in my lecture notes defined as $w_{i+1} = w_{i}+hf(t_{i}+\frac{h}{2},w_{i}+\frac{h}{2}f(t_{i},w_{i}))$ — Chengdu, Oct 04 '18 at 15:19
Ok, so the terms "improved" and "modified", as vague as they are, are probably not used consistently among the two methods. Yes, that is the midpoint method. What exactly is the problem implementing it? Following the usual interpretation of the Butcher tableau it would be `k1 = h*f(x,y); k2=h*f(x+0.5*h, y+0.5*k1); x_next=x+h; ynext = y+k2;` — Lutz Lehmann, Oct 13 '18 at 12:33

score 0 · Accepted Answer · answered Oct 04 '18 at 20:04

Below is an implementation in MATLAB I have done of the Euler's Method for solving a pair of coupled 1st order DE's. It solves a harmonic oscillator of represented by the following:

y1(t+h) = y1(t) + h*y2(t)

y2(t+h) = y2(t) + h*(-A/M y1(t) -B/M y1(t)/|y1(t)|)

% Do the integration using the Euler Method
    while(T<=T1)
        % Update the position of the pt mass using current velocity
        Y1(i+1) = Y1(i) + H*Y2(i);

        % Update the velocity of the pt mass checking if we are turning 
        % within the C/A band
        if ( (abs(Y2(i) < th) && abs(Y1(i)) < C/A) )
            Y2(i+1) = Y2(i);
        else
            Y2(i+1) = Y2(i) + H * ( ((-A/M)*Y1(i)) - (B/M)*sign(Y2(i)) );
        end

        % Incriment the time by H 
        T = T + H;
        % Increase the looping index variable
        i = i + 1;
    end

Without explicitly solving YOUR homework question, I hope this example helps. A few things to note: the if statement is accounting for static friction in this particular example -- so disregard that and only look at what happens after the 'else'.

Also note that I have initial conditions y1(0) and y2(0) defined as Y1(i=1) and Y2(i=1), so starting with Yj(i+1) gives Yj(2). Note that T1 is the ending time of the simulation.

Below is the same problem using the Improved Euler Method. If you derive the update equations for this system you should be able to walk through the code.

% Do the integration using the Improved Euler Method
    % Ki_j = K^i_j
    while(T<=T1)
        % Calculate K^i_j's

        K1_1 = Y2(i);

        % Must check if we are turning within C/A
        if ( (abs(Y2(i) < th) && abs(Y1(i)) < C/A) )
            K1_2 = 0;
        else
            K1_2 = (-A/M)*Y1(i) - (B/M)*sign(Y2(i));
        end

        K2_1 = Y2(i)+H*K1_2;

        % Checking if we are turning within C/A
        if ( (abs(Y2(i) < th) && abs(Y1(i)) < C/A) )
            K2_2 = 0;
        else
            K2_2 = (-A/M)*(Y1(i) + H*K1_1) - (B/M)*sign(Y2(i)+ H*K1_2);
        end


        % Update the position and velocity
        Y1(i+1) = Y1(i)+ (H/2)*(K1_1+K2_1);
        Y2(i+1) = Y2(i) + (H/2)*(K1_2+K2_2);

        % Incriment the time by H 
        T = T + H;
        % Increase the looping index variable
        i = i + 1;
    end

midpoint rule for matlab

1 Answers1