How to solve this?
显示 更早的评论
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?
采纳的回答
Abraham Boayue
2018-3-8
编辑：Abraham Boayue
2018-4-10
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 个评论
Abraham Boayue
2018-3-8
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
2018-3-8
Abraham Boayue
2018-3-8
You are welcome.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
