How to solve this?
6 views (last 30 days)
Show older comments
James Marlom
on 6 Mar 2018
Edited: Abraham Boayue
on 10 Apr 2018
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/171015/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/171016/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/171017/image.jpeg)
1 Comment
James Tursa
on 6 Mar 2018
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?
Accepted Answer
Abraham Boayue
on 8 Mar 2018
Edited: Abraham Boayue
on 10 Apr 2018
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
3 Comments
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!