如何用龙格库塔法来编写如下微分方程x=ax-bx^3+Acos2πft+Dy(t),a=1.b=1,f=0.01,A=0.1,D=0.5。matlab实现

Question

如题所述

Answer 1

% main.m
clc,clear,close all
t0 = 2;t1=200; h=0.02;x0=1;
T = t0:h:t1;X = zeros(length(x0),length(T));X(:,1)=x0;
for jj = 1:length(T) -1

X(:,jj+1) =RK45(T(jj),X(:,jj),@(t,x) f(t,x),h);
end
%%
figure,subplot(211),plot(T(1:100),X(1,1:100)),grid on,title('进入振子前时域')
xq = X(1,1:100);
xh = X(1,101:end);
%% 
N = length(xh);
fs = 1/h;
f = (0:N-1)*fs/N;
y = abs(fft(xh));
subplot(212),plot(f,y),grid on,title('进入振子前频域')


%%
fs =1/h;
N=length(xq);
n = 0:N-1;
y= fftshift(abs(fft(xq)));
f = n*fs/N -fs/2;

fs =1/h;
N=length(xh);
n = 0:N-1;
y= fftshift(abs(fft(xh)));
f = n*fs/N -fs/2;

figure,
subplot(211),plot(T(101:end),X(1,101:end)),grid on,title('进入振子后时域')
subplot(212),plot(f,y),grid on
title('进入振子后频域')

不知道对不对，参考下

Answer 2

x=ax-bx^3+Acos2πft+Dy(t)
式中x是自变量吗，t是什么？追问

x也是t的函数 可以写为x（t），t是时间吧，也是未知的

追答

感觉x(t)为常数处理，因为，你这里的导数为dy/dt.不知道这样理解对吗？

追问

如果x（t）设置为常数或者有初值的话又该怎么处理呢，题目中没有提到存在dy/dt的导数

本回答被网友采纳