How to solve this?

[James Tursa]

What have you done so far? What specific problems are you having with your code? Are you supposed to solve this with your own hand-written Euler and RK code, or can you use the MATLAB supplied functions such as ode45?

Answer 1

[Abraham Boayue]

function [t,x,y,N] = Runge4_2eqs(f1,f2,to,tfinal,xo,yo,h)

      N = ceil((tfinal-to)/h);
         x = zeros(1,N); 
         y = zeros(1,N) ;
         x(1) = xo;
         y(1) = yo;
         t = to;

         for i = 1:N
            t(i+1) = t(i)+h;
            Sx1 = f1(t(i),x(i),y(i));
            Sy1 = f2(t(i),x(i),y(i));

            Sx2 = f1(t(i)+h/2, x(i)+Sx1*h/2, y(i)+Sy1*h/2);
            Sy2 = f2(t(i)+h/2, x(i)+Sx1*h/2, y(i)+Sy1*h/2);

            Sx3 = f1(t(i)+h/2, x(i)+Sx2*h/2, y(i)+Sy2*h/2);
            Sy3 = f2(t(i)+h/2, x(i)+Sx2*h/2, y(i)+Sy2*h/2);

            Sx4 = f1(t(i)+h, x(i)+Sx3*h, y(i)+Sy3*h);
            Sy4 = f2(t(i)+h, x(i)+Sx3*h, y(i)+Sy3*h);

           x(i+1) = x(i) + (h/6)*(Sx1+2*Sx2+2*Sx3+Sx4);
           y(i+1) = y(i) + (h/6)*(Sy1+2*Sy2+2*Sy3+Sy4);

         end
  function f1 = DvDX(t,x,y)

% Change t to V, x to X and y to T for your problem
% To = 1035;
% vo = 0.002;
% C_Ao = 1;
% k = 3.5*exp(34222*(1/To-1/T));
% F_Ao = C_Ao*vo;
% ra = k.*C_Ao*(1-To*X)./(1+X*T);  
% f1 =  ra/F_Ao;

% Test function
% a = 10;b = 1; 
% f1= a*x-b*x*y;

f1 = y;

end
function f2 = DvDT(t,x,y) 
% Change t to V, x to X and y to T for your problem
% To = 1035;
% T_a = 1150;
% vo = .002;
% U = 110;
% a1 = 150;
% C_Ao = 1;
% F_Ao = C_Ao*vo;
% deltaHR = 80770+6.8*(T-298)-5.75e-3*(T.^2-298^2)-1.27e-6*(T.^3-298^3);
% 
% C_pA =  26.63 + 0.1830*T - (45.86e-6)*T.^2;
% C_pB =  20.04 + 0.0945*T - (30.95e-6)*T.^2;
% C_pC =  13.39 +  0.077*T - (1871e-6)*T.^2;
% 
% deltaC_p = C_pB + C_pC - C_pA;
% 
% k = 3.5*exp(34222*(1/To-1/T));
% ra = -k.*C_Ao*(1-To*X)/(1+X.*T);  
% f2 =  U*a1*(T_a-T)+ra.*deltaHR./(F_Ao*(C_pA+X.*deltaC_p));

% Test functions dydt 
% l =.1; k =1;
% f2 =  -l*y+k*x*y;

c = 0.16; m = 0.5; g = 9.81; L = 1.2;
f2 = -(c/m)*y - (g/L)*sin(x);

end
clear variables
colse all
xo = pi/2;
yo = 0;
h = .020;
to = 0;
tfinal = 20;
[t,x,y,N] = Runge4_2eqs(@DvDX,@DvDT,to,tfinal,xo,yo,h);
figure(1); clf(1) 
plot(t,x, 'Linewidth', 1.5, 'color', 'r')
hold on
plot(t,y,'Linewidth', 1.5, 'color', 'b')
legend('Dfx','Dfy')
title('Solution to two systems of ODEs')
xlabel('x')
ylabel('y')
xlim([to tfinal])
grid

[Abraham Boayue]

Hi James I posted the code above to help you get started with your exercise, it is not a solution to your problem, but a guild to help you get through it. The Runge-Kutta function is straightforward and should be easy to understand, however, due to the complexity of your equations, you might have to be very careful on selecting the step size as well as the proper selection of parameters. Euler's method can be implemented in similar fashion. Good luck for now.

[James Marlom]

Thank you Abraham

[Abraham Boayue]

You are welcome.