How to solve this?

Question

James Marlom 2018-3-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/386695-how-to-solve-this

编辑： Abraham Boayue 2018-4-10

采纳的回答： Abraham Boayue

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

James Tursa 2018-3-6

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 2018-3-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/386695-how-to-solve-this#answer_308949

编辑：Abraham Boayue 2018-4-10

在 MATLAB Online 中打开

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