清空筛选条件
清空筛选条件

Numerical differential equation solver for initial conditions not at y0

1 次查看（过去 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 中查找有关 Symbolic Math Toolbox 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by