ode45或四阶Runge-Kutta方法中的Matlab错误，用于解决耦合ODE的系统 [英] Matlab error in ode45 or fourth-order Runge-Kutta method to solve a system of coupled ODEs

问题描述

我是Matlab编程的初学者，也是Runge-Kutta方法的初学者.

I am a beginner at Matlab programming and with the Runge-Kutta method as well.

我正在尝试使用四阶Runge-Kutta方法求解耦合ODE的系统，用于我的项目工作.

I'm trying to solve a system of coupled ODEs using a 4th-order Runge-Kutta method for my project work.

这是我的问题...

G = 1.4;
g = 1.4;
k = 0;
z = 0;
b = 0.166667;

syms n; 
x2 = symfun(sym('x2(n)'),[n]);
x1 = symfun(sym('x1(n)'),[n]);
x3 = symfun(sym('x3(n)'),[n]);
x4 = symfun(sym('x4(n)'),[n]);
x5 = symfun(sym('x5(n)'),[n]);

k1 = [x2 * x1 *n *(1 - z * x2)*(x1 - n) - 2 * x3 * n *(1 - z * x2) - x4^2 * x2 *(1 - z * x2)- G *x3 *x2 ]./ [( G * x3 - (x1 - n)^2 * x2 *(1 - z * x2)) * n];   
k2 = [x2 * (1 - z * x2)*(x1 * x2 * ( x1 - 2 *n)*( x1 - n) + 2* x3 * n + x4^2 * x2 ) ]./ [( G * x3 - (x1 - n)^2 * x2 *(1 - z * x2)) * n * (x1 - n)];
k3 = [x3 * x2 * (2 * n * x1 - n)^2 * ( 1 - z * x2) + G * x1 * (x1 - 2 *n)* (x1 - n) + x4^2 * G]./ [( G * x3 - (x1 - n)^2 * x2 *(1 - z * x2)) * n * (x1 - n)];
k4 = [x4 * ( x1 + n)] ./ [n * (x1- n)];
k5 = - [x5] ./ [n * (x1- n)];

f = @(n,x) [k1;  k2;  k3;  k4; k5];
[n,xa] = ode45(f,[0 1],[1-b 1/b 1-b 0.01 0.02]);

错误是

使用odearguments时出错(第93行)
@(N，X)[K1; K2; K3; K4; K5]返回长度为1的向量，但长度为
的初始条件向量为5. @(N，X)[K1; K2; K3; K4; K5]和初始条件向量必须具有
元素数量相同.

Error using odearguments (line 93)
@(N,X)[K1;K2;K3;K4;K5] returns a vector of length 1, but the length
of initial conditions vector is 5. The vector returned by @(N,X)[K1;K2;K3;K4;K5] and the initial conditions vector must have
the same number of elements.

ode45中的错误(第114行)
[neq，tspan，ntspan，next，t0，tfinal，tdir，y0，f0，odeArgs，
odeFcn，...

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs,
odeFcn, ...

请指导我如何使用四阶Runge-Kutta方法解决上述问题...

Please guide me how can I solve the above problem with fourth-order Runge-Kutta method...

推荐答案

通常，首先请查阅文档.如果不可用，请使用您选择的搜索引擎，例如 https://www .google.de/search?q = matlab + ode45 + example 来查找第一个结果 https://de.mathworks.com/help/matlab/ref/ode45.html

In general, you consult first the documentation. If not available, use the search engine of you choice, for instance https://www.google.de/search?q=matlab+ode45+example to find as first result https://de.mathworks.com/help/matlab/ref/ode45.html

在那里您可以找到一个多维示例

There you find a multidimensional example

function dy = rigid(t,y)
    dy = zeros(3,1);    % a column vector
    dy(1) = y(2) * y(3);
    dy(2) = -y(1) * y(3);
    dy(3) = -0.51 * y(1) * y(2);
end

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T,Y] = ode45(@rigid,[0 12],[0 1 1],options);

plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')

可以用作蓝图.在相关页面中 https://de.mathworks.com/help /matlab/math/ordinary-differential-equations.html 在第一个示例中，它们对van der Pol方程使用的语法与您更接近.

that could be used as a blueprint. In the related page https://de.mathworks.com/help/matlab/math/ordinary-differential-equations.html they use a different syntax closer to yours for the van der Pol Equation in the first example

function dydt = vdp1(t,y)
    dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
end

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

plot(t,y(:,1),'-',t,y(:,2),'--')
title('Solution of van der Pol Equation, \mu = 1');
xlabel('time t');
ylabel('solution y');
legend('y_1','y_2')

按照以下示例，将您的代码重写为

Following these examples, rewrite your code as

G = 1.4;
g = 1.4;
k = 0;
z = 0;
b = 0.166667;



function dotx = dxdn(n,x) 
    t1 = n*(x1-n)    
    t2 = x(2)*(1 - z * x(2))

    den123 = ( G * x(3) - (x(1) - n)^2 * t2)

    k1 = ( x(1)*t1*t2 - 2 * x(3) * n *(1 - z * x(2)) - x(4)^2 * t2- G *x(3) *x(2) ) / (den123*n);   
    k2 = ( t2*(x(1) * x(2) * ( x(1) - 2 *n)*( x(1) - n) + 2* x(3) * n + x(4)^2 * x(2) ) ) / (den123 * t1);
    k3 = ( x(3) *  (2 * n * x(1) - n)^2 * t2 + G * x(1) * (x(1) - 2 *n)* (x(1) - n) + x(4)^2 * G ) / (den123*t1);
    k4 = ( x(4) * ( x(1) + n) ) / t1;
    k5 = - x(5) / t1;

    dotx = [k1;  k2;  k3;  k4; k5];
end 

[n,xa] = ode45(@dxdn,[0.001 1],[1-b; 1/b; 1-b; 0.01; 0.02]);

由于在n=0处被零除，因此无法以这种奇异性开始迭代.在上面的代码中，可以通过以(非常)小的正数n开头来缓解这种情况，也可以尝试以n=1e-8或更小数开头.所有组件的斜率都非常大，因此积分可能会很慢，并且结果可能不会过于精确地接近零.要正确处理奇异的ODE，请在math.stackexchange论坛中提出.

Because of the division by zero at n=0 you can not start the iteration in that singularity. In the code above, this is mitigated by starting with some (very) small positive n, you can also try starting with n=1e-8 or smaller. The slope will be very large in all components, so the integration may be slow, and the result might be not overly exact close to zero. For the correct handling of singular ODE ask in the math.stackexchange forum.

这篇关于ode45或四阶Runge-Kutta方法中的Matlab错误，用于解决耦合ODE的系统的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！