Solving system of ODEs with ode45

조회 수: 1 (최근 30일)
shir levi
shir levi 2022년 12월 24일
댓글: Sam Chak 2022년 12월 25일
i have the following equetions system
and and the maching valuse
i tried to solve the system with the following code but i dont get the right resultes which i know:
ephsilon values [1.94 2.65 3.42] for the x1/h vec [7.5 9.5 11]
k values [0.2680 0.3585 0.444] for the same x1/h vec.
the code i wrote
Y=zeros(1,2);
Uc=12.4; %m/s
S=46.8;% 1/s
k_75=0.268 ; %m2/s2;
e_75=1.94; %m2/s3
Cu = 0.09;
C1= 1.44;
C2=1.92;
x1_h=7.5;
h = 0.3048;%m
xspan=[7.5 9.5 11];
Y0=[k_75 e_75 ];
f = @(t,Y) [(Cu*((Y(1))^2)/(Y(2))*(S^2)-Y(2))/(Uc/h); (Cu*C1*(Y(1))*(S^2)-C2*((Y(2))^2)/(Y(1)))/(Uc/h)];
[x0,val] = ode45(f , xspan, Y0);
figure (1)
plot(x0, val(:,1))
figure (2)
plot(x0, val(:,2))
what wrong with it ? thanks for the help
  댓글 수: 1
Walter Roberson
Walter Roberson 2022년 12월 25일
What is the evidence that those are not correct output graphs?

댓글을 달려면 로그인하십시오.

답변 (2개)

Sam Chak
Sam Chak 2022년 12월 24일
편집: Sam Chak 2022년 12월 24일
Edit: Previously I've got the same results like you did. However, after thinking about it, those 2 sets of 3 values for k and ϵ are probably initial conditions. So, I modified the code to include a for loop.
Also, in the 3rd figure, the system is shown to be inherently unstable, and there is no equilibrium point.
% parameters
Uc = 12.4; % m/s
S = 46.8; % 1/s
k = [0.2680 0.3585 0.444]; % m2/s2 ic for k
epsil = [1.9400 2.6500 3.420]; % m2/s3 ic for epsil
Cu = 0.09;
C1 = 1.44;
C2 = 1.92;
x1_h = [7.5 9.5 11];
h = 0.3048; % m
tstart = x1_h*h/Uc;
tfinal = 2*tstart;
for j = 1:length(x1_h)
% setting conditions
f = @(t, Y) [Cu*(Y(1)^2)/Y(2)*S^2 - Y(2); Cu*C1*Y(1)*S^2 - C2*(Y(2)^2)/Y(1)];
tspan = [tstart(j) tfinal(j)];
Y0 = [k(j) epsil(j)];
[t, Y] = ode45(f, tspan, Y0);
% plotting results
figure(1)
plot(t*Uc/h, Y(:,1)), hold on, grid on,
xlabel('x_{1}/h'), ylabel('k'), title('k responses')
figure(2)
plot(t*Uc/h, Y(:,2)), hold on, grid on,
xlabel('x_{1}/h'), ylabel('\epsilon'), title('\epsilon responses')
end
hold off
[x, y] = meshgrid(-11:2:11);
u = Cu*(x.^2)./y*S^2 - y;
v = Cu*C1*x*S^2 - C2*(y.^2)./x;
figure(3)
l = streamslice(x, y, u, v);
axis tight
xlabel('k')
ylabel('\epsilon')

shir levi
shir levi 2022년 12월 24일
thak you for answaring
but those arnt the initial conditions
its the valuse i was supposed to get in my results
  댓글 수: 1
Sam Chak
Sam Chak 2022년 12월 25일
But you used these values and in the initial condition.
So, I'm confused.
Perhaps, if you tell us more info about the 2 sets of 3 values, it will help clarifying the matter.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by