Numerical differential equation solver for initial conditions not at y0

2 次查看（过去 30 天）
Rihards
Rihards 2022-11-6
评论： Torsten 2022-11-6
I've been trying to convert a part of mathematica code to a Matlab one, in order to integrate it for other calculations, but I am having some trouble getting the solver to run.
If by any chance this helps, this is the script that I am trying to convert:
And this is my approach.
syms Uf(t) delta y
d = 0.2000;
k_Karman = 0.4;
sol_Um = 3.2947;
omega = 0.5405;
Tp = 29.0823;
ks = 2.5 * d *10^-3 % [m]
k_Karman = 0.4 % von Karman constant
U_R04 = sol_Um * sin(omega*t)
U0_R04 = Uf/k_Karman*log(30*y/ks)
eq_R04 = sol_Um * sin(omega*t) == Uf/k_Karman*log(30*delta/ks)
sol_delta = solve(eq_R04,delta) % solve for delta
ME_R04 = int(diff(U0_R04,t),y,[0 sol_delta]) == sol_delta*diff(U_R04,t)-Uf^2
I believe that this far, I would be getting a matching result from the example below, however, I am unable to prepare the command line for the evaluation.
This was my initial attempt, but I have a feeling that is the wrong approach to take for this specific case as I am getting an error for the 'Vars' definition.
[VF,Sbs] = odeToVectorField(ME_R04);
Uffcn = matlabFunction(VF, 'Vars',{x,Y});
[x,Ufx] = ode15s(Uffcn, [0 Tp/2], 0);
Error using sym/matlabFunction>getOptions
The value of 'Vars' is invalid. 'Vars' value must be a character vector, a 1-dimensional cell array of character vectors, a 1-dimensional cell array of symbolic variables
or arrays of symbolic variables, or an array of symbolic variables.
Initial attempt was done for initial conditions of y0=0, as I am not aware of a way to define it the way mathematica script allows for Uf(0.001*Tp)=0.001*sol_Um
  1 个评论
Rihards
Rihards 2022-11-6
Edited the last part as follows:
[VF,Sbs] = odeToVectorField(ME_R04);
Uffcn = matlabFunction(VF, 'Vars',{t,Y});
init_cond = 0.0032947071555424111234342753369588 % Gives error when defined as 0.001*Tp
[t,Ufx] = ode15s(Uffcn, [0.001*Tp Tp/2],init_cond);
figure(9);
plot(t, Ufx,'b','linewidth',1)
And I am curious if anyone can explain, why it solves here, but gives an error when the init_cond is defined as
init_cond = 0.001*Tp

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-11-6
Ran in:
The solution seems to be very sensitive to the initial condition.
syms Uf(t) delta y
d = 0.2000;
k_Karman = 0.4;
sol_Um = 3.2947;
omega = 0.5405;
Tp = 29.0823;
ks = 2.5 * d *10^-3; % [m]
k_Karman = 0.4; % von Karman constant
U_R04 = sol_Um * sin(omega*t);
U0_R04 = Uf/k_Karman*log(30*y/ks);
eq_R04 = sol_Um * sin(omega*t) == Uf/k_Karman*log(30*delta/ks);
sol_delta = solve(eq_R04,delta); % solve for delta
Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.
ME_R04 = int(diff(U0_R04,t),y,[0 sol_delta]) == sol_delta*diff(U_R04,t)-Uf^2;
[VF,Sbs] = odeToVectorField(ME_R04);
Uffcn = matlabFunction(VF);
init_cond = 0.02;%0.001*Tp%0.0032947071555424111234342753369588; % Gives error when defined as 0.001*Tp
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
[t,Ufx] = ode15s(@(t,y)Uffcn(y,t), [0.001*Tp 5.8],init_cond,options);
figure(1);
plot(t, Ufx,'b','linewidth',1)
  2 个评论
Rihards
Rihards 2022-11-6
I am getting very similar result with this. Is there a way to increase the number of solutions for the plot? Currently I got 72 values for t and Ufx, but I believe the very peak value is somewhere in between them.
syms Uf(t) delta y
ks = 2.5 * d *10^-3 % [m]
k_Karman = 0.4 % von Karman constant
U_R04 = sol_Um * sin(omega*t)
U0_R04 = Uf/k_Karman*log(30*y/ks)
eq_R04 = sol_Um * sin(omega*t) == Uf/k_Karman*log(30*delta/ks)
sol_delta = solve(eq_R04,delta) % solve for delta
ME_R04 = int(diff(U0_R04,t),y,[0 sol_delta]) == sol_delta*diff(U_R04,t)-Uf^2
[VF,Sbs] = odeToVectorField(ME_R04);
Uffcn = matlabFunction(VF, 'Vars',{t,Y});
init_cond = 0.0032947071555424111234342753369588 % Gives error when defined directly as 0.001*sol_Um
[t,Ufx] = ode15s(Uffcn, [0.001*T T/2],init_cond);
figure(9);
plot(t, Ufx,'b','linewidth',1)
title('R04-Task2','interpreter','latex','fontsize',16)
xlabel('Time, [s]','interpreter','latex','fontsize',16)
ylabel('Friction velocity Uf, [m/s]','interpreter','latex','fontsize',16)
grid on
Torsten
Torsten 2022-11-6
nt = 150;
[t,Ufx] = ode15s(Uffcn, linspace(0.001*T T/2,nt),init_cond);
instead of
[t,Ufx] = ode15s(Uffcn, [0.001*T T/2],init_cond);

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by