Matlab Numerical Solver:求解二阶微分方程
Matlab Numerical Solver: Solving Second Order differential Equation
非常抱歉,我无法提供我正在使用的确切方程式的详细信息。这是一个非常复杂的二阶微分方程,形式类似于:
其中函数 a(z) ~ e(z) and g(z)
给出。 p
不变。
边界条件我也有
是否可以借助matlab求解f(z)
?
如有任何建议,我们将不胜感激。非常感谢。
编辑
这是我定义函数的方式。 a1 ~ g1
和 fna ~ fng
定义并存储在 Gdata.mat
:
function xp=myfunc(t,p)
% x = [d2f df f df2 f2 d2f2]
% xp = [df d2f f2 df2 d3f d2f2]
load GData
xp=zeros(6,1); % [f df d2f f2 df2]
% f
fprintf('%d\n',length(xp));
fprintf('%d\n',length(p));
xp(1) = x(2); % df
xp(2) = x(1); % d2f
% f2
xp(4) = x(4); % df2
xp(6) = x(6); % d2f2
xp(5) = (...
b1(t)*p(3) + b(t)*p(2) + ...
c1(t)*p(3)^3 + 3*fnc(t)*p(3)^2*p(2) + ...
d1(t)*p(3)^5 + 5*fnd(t)*p(3)^4*p(2) + ...
e1(t)*p(3)^7 + 7*fne(t)*x(3)^6*p(2) - ...
f1(t)*p(2)*p(3) + f1(t)*p(1)*p(3) + f1(t)*p(2)^2 - ...
g1(t)*p(4) - fng(t)*p(6) + ...
q*p(2) - a1(t)*p(1)...
) * 1/(fna(t));
然后我打电话给:
[TEMP,POL] = ode45('odesolver',[0,1],[0,0,0,0,0,0]);
编辑2
function dp=odesolver(t,p)
% dp = [df d2f d3f]
syms x;
load BData;
load GData;
dp=zeros(3,1); % [f df d2f]
A = interp1(t_data,At,t);
B = interp1(t_data,Bt,t);
C = interp1(t_data,Ct,t);
D = interp1(t_data,Dt,t);
E = interp1(t_data,Et,t);
F = interp1(t_data,Ft,t);
G = interp1(t_data,Gt,t);
A1 = interp1(t_data,A1t,t);
B1 = interp1(t_data,B1t,t);
C1 = interp1(t_data,C1t,t);
D1 = interp1(t_data,D1t,t);
E1 = interp1(t_data,E1t,t);
F1 = interp1(t_data,F1t,t);
G1 = interp1(t_data,G1t,t);
dp(1) = p(2); % f'
dp(2) = p(3); % f''
dp(3) = (...
B1*p(3) + B*p(2) + ...
C1*p(3)^3 + 3*C*p(3)^2*p(2) + ...
D1*p(3)^5 + 5*D*p(3)^4*p(2) + ...
E1*p(3)^7 + 7*E*p(3)^6*p(2) - ...
F1*p(2)*p(3) + F*p(1)*p(3) + F*p(2)^2 - ...
2*G1*p(2)*p(1) + 2*G*p(3)*p(1) + 2*G*p(2)*p(2) + ...
q*p(2) - A1*p(1)) * 1/(A);
根据讨论和编辑过的问题回答:
使用 ode45
求解微分方程有几个障碍,但 none 是个阻碍:
- 二阶 ODE:转换为 2 个一阶方程,您可以使用
ode45
求解,如 this question.
- 积分项:微分你的方程来摆脱它。您最终会得到一个三阶微分方程,您需要使用与上述相同的技术将其转换为 3 个一阶方程
- 要积分的函数的平方的导数:向您的
ode
函数添加额外的状态以允许使用ode45
求解
- 时间相关变量:你不能直接使用你的时间变量来索引你的时间相关变量。您需要使用随时间变化的数据对
t
的每个值进行线性插值以获得变量的相应值,如 this other question.
所以你的代码应该看起来像这样(我假设你的 GData
文件中的所有时间相关数据都有一个共同的时间向量,这可能不是这种情况,你可能有每个变量的不同时间向量 - 如有必要调整代码;无论哪种方式,请确保您的时间相关数据是在您调用的时间间隔上定义的 ode45
):
function dx=myfunc(t,x)
% x = [d2f df f df2 f2 d2f2]
% dx = [d3f d2f df d2f2 df2 d3f2]
load GData % I assume this contains b1, b, c1, etc... *and* the corresponding time vector, say t_data
dx=zeros(6,1); % [d3f d2f df d2f2 df2 d3f2]
fprintf('%d\n',length(dx));
fprintf('%d\n',length(x));
dx(1) = x(2); % double-check!!
dx(2) = x(1); % double-check!!
% double-check!!
dx(4) = x(4); % double-check!!
dx(6) = x(6); % double-check!!
B1 = interp1(t_data,b1,t);
B = interp1(t_data,b,t);
C1 = interp1(t_data,c1,t);
% etc... for the other variables
dx(5) = (...
B1*p(3) + B*p(2) + ...
C1*p(3)^3 + 3*FNC*p(3)^2*p(2) + ...
D1*p(3)^5 + 5*FND*p(3)^4*p(2) + ...
E1*p(3)^7 + 7*FNE*x(3)^6*p(2) - ...
F1*p(2)*p(3) + F1*p(1)*p(3) + F1*p(2)^2 - ...
G1*p(4) - FNG*p(6) + ...
q*p(2) - A1*p(1)...
) * 1/(FNA); % double-check!!
非常抱歉,我无法提供我正在使用的确切方程式的详细信息。这是一个非常复杂的二阶微分方程,形式类似于:
a(z) ~ e(z) and g(z)
给出。 p
不变。
边界条件我也有
是否可以借助matlab求解f(z)
?
如有任何建议,我们将不胜感激。非常感谢。
编辑
这是我定义函数的方式。 a1 ~ g1
和 fna ~ fng
定义并存储在 Gdata.mat
:
function xp=myfunc(t,p)
% x = [d2f df f df2 f2 d2f2]
% xp = [df d2f f2 df2 d3f d2f2]
load GData
xp=zeros(6,1); % [f df d2f f2 df2]
% f
fprintf('%d\n',length(xp));
fprintf('%d\n',length(p));
xp(1) = x(2); % df
xp(2) = x(1); % d2f
% f2
xp(4) = x(4); % df2
xp(6) = x(6); % d2f2
xp(5) = (...
b1(t)*p(3) + b(t)*p(2) + ...
c1(t)*p(3)^3 + 3*fnc(t)*p(3)^2*p(2) + ...
d1(t)*p(3)^5 + 5*fnd(t)*p(3)^4*p(2) + ...
e1(t)*p(3)^7 + 7*fne(t)*x(3)^6*p(2) - ...
f1(t)*p(2)*p(3) + f1(t)*p(1)*p(3) + f1(t)*p(2)^2 - ...
g1(t)*p(4) - fng(t)*p(6) + ...
q*p(2) - a1(t)*p(1)...
) * 1/(fna(t));
然后我打电话给:
[TEMP,POL] = ode45('odesolver',[0,1],[0,0,0,0,0,0]);
编辑2
function dp=odesolver(t,p)
% dp = [df d2f d3f]
syms x;
load BData;
load GData;
dp=zeros(3,1); % [f df d2f]
A = interp1(t_data,At,t);
B = interp1(t_data,Bt,t);
C = interp1(t_data,Ct,t);
D = interp1(t_data,Dt,t);
E = interp1(t_data,Et,t);
F = interp1(t_data,Ft,t);
G = interp1(t_data,Gt,t);
A1 = interp1(t_data,A1t,t);
B1 = interp1(t_data,B1t,t);
C1 = interp1(t_data,C1t,t);
D1 = interp1(t_data,D1t,t);
E1 = interp1(t_data,E1t,t);
F1 = interp1(t_data,F1t,t);
G1 = interp1(t_data,G1t,t);
dp(1) = p(2); % f'
dp(2) = p(3); % f''
dp(3) = (...
B1*p(3) + B*p(2) + ...
C1*p(3)^3 + 3*C*p(3)^2*p(2) + ...
D1*p(3)^5 + 5*D*p(3)^4*p(2) + ...
E1*p(3)^7 + 7*E*p(3)^6*p(2) - ...
F1*p(2)*p(3) + F*p(1)*p(3) + F*p(2)^2 - ...
2*G1*p(2)*p(1) + 2*G*p(3)*p(1) + 2*G*p(2)*p(2) + ...
q*p(2) - A1*p(1)) * 1/(A);
根据讨论和编辑过的问题回答:
使用 ode45
求解微分方程有几个障碍,但 none 是个阻碍:
- 二阶 ODE:转换为 2 个一阶方程,您可以使用
ode45
求解,如 this question. - 积分项:微分你的方程来摆脱它。您最终会得到一个三阶微分方程,您需要使用与上述相同的技术将其转换为 3 个一阶方程
- 要积分的函数的平方的导数:向您的
ode
函数添加额外的状态以允许使用ode45
求解
- 时间相关变量:你不能直接使用你的时间变量来索引你的时间相关变量。您需要使用随时间变化的数据对
t
的每个值进行线性插值以获得变量的相应值,如 this other question.
所以你的代码应该看起来像这样(我假设你的 GData
文件中的所有时间相关数据都有一个共同的时间向量,这可能不是这种情况,你可能有每个变量的不同时间向量 - 如有必要调整代码;无论哪种方式,请确保您的时间相关数据是在您调用的时间间隔上定义的 ode45
):
function dx=myfunc(t,x)
% x = [d2f df f df2 f2 d2f2]
% dx = [d3f d2f df d2f2 df2 d3f2]
load GData % I assume this contains b1, b, c1, etc... *and* the corresponding time vector, say t_data
dx=zeros(6,1); % [d3f d2f df d2f2 df2 d3f2]
fprintf('%d\n',length(dx));
fprintf('%d\n',length(x));
dx(1) = x(2); % double-check!!
dx(2) = x(1); % double-check!!
% double-check!!
dx(4) = x(4); % double-check!!
dx(6) = x(6); % double-check!!
B1 = interp1(t_data,b1,t);
B = interp1(t_data,b,t);
C1 = interp1(t_data,c1,t);
% etc... for the other variables
dx(5) = (...
B1*p(3) + B*p(2) + ...
C1*p(3)^3 + 3*FNC*p(3)^2*p(2) + ...
D1*p(3)^5 + 5*FND*p(3)^4*p(2) + ...
E1*p(3)^7 + 7*FNE*x(3)^6*p(2) - ...
F1*p(2)*p(3) + F1*p(1)*p(3) + F1*p(2)^2 - ...
G1*p(4) - FNG*p(6) + ...
q*p(2) - A1*p(1)...
) * 1/(FNA); % double-check!!