MATLAB的答案

拟合二阶微分方程参数

25次浏览(过去30天)
克里斯托弗·奈尔
克里斯托弗·奈尔 2021年9月22日
评论道: 明星黾 2021年9月24日
我有一个带有可调参数的二阶微分方程我需要用它来拟合一些数据。我花了一些时间看了一些其他帖子的答案(特别是StarStrider, Torsten的一些关键助攻)，最终让我的脚本运行，但有些地方不正确。我怀疑这是我在将y(t)数据传递到包含微分方程组的函数时遗漏的东西，但我不确定。我希望有人能看到我错过了什么。数据拟合函数应该更紧密地跟随数据。数据是y vs t，而div . eq是y''=a，其中a是y'的函数。以下是我的剧本:
清晰的所有
关闭所有
clc
Time_data = [0.8273;1.1104;1.3940;1.6584;1.8764;2.0173;2.1239;2.2343;2.3498) -。35;%自变量，实验数据
Y_data = [1.5025;2.4650;3.4274;4.3976;5.3600;6.3302;7.2965;8.2590;9.2253];%因变量，实验数据
A_m_guess = 200;a_m的%init猜测，要拟合的参数
a_m = lsqcurvefit(@(am,t)FittingFunction(am,t)，a_m_guess,Time_data,Y_data);
数字绘制数据以及拟合函数
持有“上”
情节(Time_data Y_data,“柯”，“DisplayName的”，“数据”）
plot_times = linspace(0,max(Time_data)，100);
情节(plot_times FittingFunction (a_m plot_times),“- r”，“DisplayName的”，“健康”);
包含(“时间”), ylabel (“Y”),传说(“位置”，“西北”）
函数F = FittingFunction(a_m,Time_data)
Init_slope = 3.39;这是一个边界条件
tspan =时间数据;
InitConds = [0 init_slope];% y (t = 0) = 0, y ' (t = 0) =
[t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,a_m)，tspan,InitConds);
F = y_atTdata(:，1);%只返回y，不返回y和y'
结束
函数Y = SysOfDiffEqns(t, Y,a_m)
% y''=a(a_m,y')是二阶差分eq。
% a是a_m()和y'(t)的函数
%需要使a_m适合可用的t_data和y_data对
使一阶deff方程组为，
%替换Y(1) = Y
% Y(2) = Y '
%生成两个一阶ode的系统:
% y (1)' = y (2)
% Y(2)' = a
Y = 0 (2,1);初始化输出矩阵
D = 8.7049;%常数输入参数
U_m = 0.95*D;%常量输入参数
U = y(2);% U等于dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m));%加速度函数，仅对U
A = max(A,0);%的力至少为零
Y(1) = Y(2);%一阶导数替换
Y(2) = a;%二阶导数等于a (y'的函数)
结束
非常感谢您的帮助。同样，我在想我需要以某种方式将我的y_data传递到@SysOfDiffEqns，然后设置y(1)=y_data，但我无法在摆弄它几天后弄清楚。

登录回答这个问题。

接受的答案

明星黾
明星黾 2021年9月22日
跑:
问题似乎是不允许初始条件作为参数以及的初始值进行估计 “a_m” ．
同时, “plot_Times” 必须仅在拟合的时间范围内定义。如果你想推断 0 (或任何其他值)，则使用 德瓦尔 函数。
我唯一做的改变就是 “a_m_guess” 在脚本和 “FittingFunction” 函数, “InitConds” 赋值和单个参数 “a_m(1)” 传递给 “SysOfDiffEqns” 在 数值 调用。
只有这些更改，代码似乎工作，并给出适当的结果-
Time_data = [0.8273;1.1104;1.3940;1.6584;1.8764;2.0173;2.1239;2.2343;2.3498) -。35;%自变量，实验数据
Y_data = [1.5025;2.4650;3.4274;4.3976;5.3600;6.3302;7.2965;8.2590;9.2253];%因变量，实验数据
A_m_guess = [200;1;1);a_m的%init猜测，要拟合的参数
a_m = lsqcurvefit(@(am,t)FittingFunction(am,t)，a_m_guess,Time_data,Y_data)
局部最小值。Lsqcurvefit停止，因为当前步长的大小小于步长公差的值。
a_m = 3×1
369.4651 1.0532 4.2399
数字绘制数据以及拟合函数
持有“上”
情节(Time_data Y_data,“柯”，“DisplayName的”，“数据”）
plot_times = linspace(min(Time_data)，max(Time_data)，100);
情节(plot_times FittingFunction (a_m plot_times),“- r”，“DisplayName的”，“健康”);
包含(“时间”), ylabel (“Y”),传说(“位置”，“西北”）
函数F = FittingFunction(a_m,Time_data)
Init_slope = 3.39;这是一个边界条件
tspan =时间数据;
% InitConds = [0 init_slope];% y (t = 0) = 0, y ' (t = 0) =
InitConds = a_m(2:3);
[t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,a_m(1))，tspan,InitConds);
F = y_atTdata(:，1);%只返回y，不返回y和y'
结束
函数Y = SysOfDiffEqns(t, Y,a_m)
% y''=a(a_m,y')是二阶差分eq。
% a是a_m()和y'(t)的函数
%需要使a_m适合可用的t_data和y_data对
使一阶deff方程组为，
%替换Y(1) = Y
% Y(2) = Y '
%生成两个一阶ode的系统:
% y (1)' = y (2)
% Y(2)' = a
Y = 0 (2,1);初始化输出矩阵
D = 8.7049;%常数输入参数
U_m = 0.95*D;%常量输入参数
U = y(2);% U等于dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m));%加速度函数，仅对U
A = max(A,0);%的力至少为零
Y(1) = Y(2);%一阶导数替换
Y(2) = a;%二阶导数等于a (y'的函数)
结束
做出适当的改变以获得不同的结果。
．
8的评论
显示旧的评论
明星黾
明星黾 2021年9月24日
跑:
我致力于运行代码并发布它，以造福于遇到这种情况并希望看到结果的任何人-
% ydot_o = 3.39;%初始斜率估计值
% D_guess = 8.7049;最终斜率的%估计
Time_data = [0.8273;1.1104;1.3940;1.6584;1.8764;2.0173;2.1239;2.2343;2.3498];%自变量，实验数据
Y_data = [1.5025;2.4650;3.4274;4.3976;5.3600;6.3302;7.2965;8.2590;9.2253];%因变量，实验数据
注意，数据不是从t=0开始的，所以初始条件是数据开始的时间
Guess = [1000, 1.5, 3.3, 8.7, 0.94];% [a_m, y(t_o)， y'(t_o)， D, U_m_factor] init guess for a_m, y_o, and y'_o, D, Um factor
Constrain_parameters = 1;
如果Constrain_parameters == 0
Lower_bound = [];
Upper_bound = [];
elseifConstrain_parameters == 1
Lower_bound = [100, 1,3,8,0.91];
Upper_bound = [1e4, 2,4,9,0.98];
结束
Return_y_and_ydot = 0;% if =0，只返回y值(用于拟合)。如果=1，返回y和y'值(用于绘图)。
Options = optimset();
[拟合，resnorm, residuals，~，~] = lsqcurvefit(@(拟合，t)FittingFunction(拟合，t,return_y_and_ydot)，guess,Time_data,Y_data,lower_bound,upper_bound,options)
局部最小值。Lsqcurvefit停止了，因为相对于其初始值的平方和的最终变化小于函数公差的值。
安装= 1×5
1.0e+03 * 1.1093 0.0015 0.0032 0.0087 0.0009
Resnorm = 0.0027
残差= 9×1
0.0131 -0.0226 -0.0019 0.0205 -0.0099 0.0204 -0.0183 -0.0201 0.0187
%%计算和绘图
Return_y_and_ydot = 1;
Trans_threshold = 0.99;
D =拟合(4);
trans_pt = trans_threshold*D;
图(1)
次要情节(2,1,1)绘制y(t)数据以及拟合的y(t)
持有“上”
情节(Time_data Y_data,“柯”，“DisplayName的”，“数据”）
plot_times = linspace(min(Time_data)，max(Time_data)，1e4);
解决方案= FittingFunction(拟合，plot_times,return_y_and_ydot);
情节(plot_times、解决方案(:1),“- r”，“DisplayName的”，“健康”);
包含(“时间”), ylabel (“Y”),传说(“位置”，“西北”)、标题(“数据和契合度vs时间”)、网格“上”
次要情节(2,1,2)将拟合的y'(t)与过渡线绘制
持有“上”
情节(plot_times、解决方案(:,2),“- r”，“DisplayName的”，“健康”);
线([min (Time_data), max (Time_data)], [trans_pt trans_pt],“DisplayName的”，“过渡”）
包含(“时间”), ylabel (“Ydot”),传说(“位置”，“西北”)、标题(“ydot vs time”)、网格“上”
图(2)%表示残差
情节(Time_data残差,“柯”，“DisplayName的”，“残差”）
包含(“时间”), ylabel (“Y”),传说(“位置”，“西北”)、标题(“剩余收入vs时间”)、网格“上”
%%向后推断解到t=0
ext_times = linspace(0,min(Time_data)，100);%延长时间到0
ICs =[拟合(2)，拟合(3)];新的集成电路是什么?
ext_sol = ode23s(@(t,y)SysOfDiffEqns(t,y, fitting)，[min(Time_data)， 0]，ICs);
Ext_y = deval(ext_sol,ext_times);
图(1)
次要情节(2,1,1)% y(t)图
情节(ext_times ext_y (1:)“g”，“DisplayName的”，“扩展”）
次要情节(2,1,2)% y'(t)图
:情节(ext_times ext_y (2),“g”，“DisplayName的”，“扩展”）
% %功能
函数F = FittingFunction(guess,Time_data,return_y_and_ydot)
tspan =时间数据;
InitConds = [guess(2)， guess(3)];% y (t = t_o), y ' (t = t_o) = init_slope
% [t_atTdata,y_atTdata] = ode45(@(t,y)SysOfDiffEqns(t,y,guess)，tspan,InitConds);
[t_atTdata,y_atTdata] = ode23s(@(t,y)SysOfDiffEqns(t,y,guess)，tspan,InitConds);
如果Return_y_and_ydot == 0%用于对数据进行优化操作
F = y_atTdata(:，1);%只返回y，不返回y和y'
elseifReturn_y_and_ydot == 1%用于显示结果
F = y_atTdata;%返回y和y'
结束
结束
函数Y = SysOfDiffEqns(t, Y,guess)
% y''=a(a_m,y')是二阶差分eq。
% a是a_m()和y'(t)的函数
%需要使a_m适合可用的t_data和y_data对
使一阶deff方程组为，
%替换Y(1) = Y
% Y(2) = Y '
%生成两个一阶ode的系统:
% y (1)' = y (2)
% Y(2)' = a
Y = 0 (2,1);初始化输出矩阵
A_m = guess(1);
D = guess(4);%让参数拟合
U_m_factor = guess(5);%让参数拟合
% d = 8.7049;%常数输入参数
% U_m_factor = 0.94;%常数输入参数
U_m = U_m_factor*D;
U = y(2);% U等于dy/dt
a = a_m*((D-U)/(D-U_m))*exp((U-U_m)/(D-U_m));%加速度函数，仅对U
A = max(A,0);%的力至少为零
Y(1) = Y(2);%一阶导数替换
Y(2) = a;%二阶导数等于a (y'的函数)
结束
我祝贺扩展原来的想法!
．

登录评论。

更多答案(0)

类别

了解更多假设在帮助中心而且文件交换

标签

s manbetx 845


释放

R2017b

社区寻宝

在MATLAB Central中找到宝藏，并发现社区如何帮助您!

开始狩猎!