我想用自己的runge kutty方法在matlab中求解以下的偏微分方程组: $/frac{dx_{1}{dt}$ = $x_{2}$ $/frac{dx_{2}{dt}$ = A $x_{1}。+ B $x_{2}$ + asin(C)t) 我写道
function xdot = NaszModel(x, aa, bb, AA, BB, CC)
a= aa;
A = AA;
B = BB;
C = CC;
xdot = [ x(2)
A*x(1) - B*x(2) + a*sin(C)];
end
function [x, t] = RK(f, x0, t0, T, h)
t = t0:h:T;
nt = numel(t);
nx = numel(x0);
x = nan(nx, nt);
x(:,1) = x0;
for k = 1:nt-1
k1 = h*f(t(k), x(:,k));
k2 = h*f(t(k) + h/3 , x(:,k) + h/3 );
k3 = h*f(t(k) + (2*h)/3, x(:,k) - k1/3 + k2);
k4 = h*f(t(k) + h, x(:,k) + k1 - k2 + k3);
dx = (k1 + 3*k2 + 3*k3 + k4)/8;
x(:,k+1) = x(:,k) + dx;
end
end
a = 1;
A = 0.1;
B = 0.1;
C = 1;
f = @(t,x) NaszModel(x, a, b, A, B, C);
x0 = [1, 0];
t0 = 0;
T = 100;
h = 1;
[x, t] = RK(f, x0, t0, T, h);
plot(t,x);
legend('x', 'v');
plot(x(1,:), x(2,:));
xlabel('x');
ylabel('v(x)');
一切都很好,但当我试图将t加入到sin中时,我的意思是:我的解决方案中包含一个sin(C),但我必须用一个公式来解决sin(C*t)。`
应该如何写才能使它也能适用于一个人sin(C)t)?
你需要通过 t
的模型函数,那么你也可以在 ODE 中使用它。这意味着你需要改变函数的接口。
function xdot = NaszModel(t, x, a, b, A, B, C)
xdot = [ x(2)
A*x(1) - B*x(2) + a*sin(C*t)];
end
f = @(t,x) NaszModel(t, x, a, b, A, B, C);