Matlab Numerical Solver:求解二阶微分方程

Matlab Numerical Solver: Solving Second Order differential Equation

非常抱歉,我无法提供我正在使用的确切方程式的详细信息。这是一个非常复杂的二阶微分方程,形式类似于:

其中函数 a(z) ~ e(z) and g(z) 给出。 p 不变。

边界条件我也有

是否可以借助matlab求解f(z)

如有任何建议,我们将不胜感激。非常感谢。

编辑

这是我定义函数的方式。 a1 ~ g1fna ~ 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!!