Solving system of ODEs with ode45

1 view (last 30 days)
shir levi
shir levi on 24 Dec 2022
Commented: Sam Chak on 25 Dec 2022
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 Comment
Walter Roberson
Walter Roberson on 25 Dec 2022
What is the evidence that those are not correct output graphs?

Sign in to comment.

Answers (2)

Sam Chak
Sam Chak on 24 Dec 2022
Edited: Sam Chak on 24 Dec 2022
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 on 24 Dec 2022
thak you for answaring
but those arnt the initial conditions
its the valuse i was supposed to get in my results
  1 Comment
Sam Chak
Sam Chak on 25 Dec 2022
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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!