matlab - Matlabのルンゲ・クッタ連続法による運動方程式

Question

• 運動方程式は次の式で与えられます。

、ここで、m、b は質量と減衰の定常値です。時変項 f(t) は励起電力で、q(t) は一般化された変位です。

• 私はそれを解決しました：

• そして、[t、x] = ode23（ 'rightside'、tspan、x0）を介してMatLabで解決する必要があります。
• f(t) と k(t) は、Matlab のようなフーリエ級数を介して複素数で解決しました。

% f(t) の複素フーリエ級数

ft=zeros(size(t)); 
for j=1:2*N+1
n= j-(N+1);    
if n==0
   f(j)=f0/2;
else
   f(j)=f0*( (exp(-i*n*2*pi)*(i*2*pi*n+1)-1)/(4*pi^2*n^2));   
   end
ft=ft+f(j)*exp(i*n*om*t);
end

% k(t) の複素フーリエ級数

kt=k0*ones(size(t)); 
for s=1:2*N+1
    n= s-(N+1);    
    if n==0
       c(s)=k0;
    else
       c(s)=i*(k0+ktyl)/n/pi*(1-cos(n*pi));     
    end
    kt=kt+c(s)*exp(i*n*om*t);
end

• そして、私たちは知っています：

T=30;
dt=0.01;
t=0:0.01:5*T;
k0=1e6;
om=2*pi/T;
ktyl=0.5e6;
N=10;
m=1;
ks=1e4;
D=0.01;
OMG=sqrt(ks/m);
b=2*D*OMG*m;
f0=100;

ありがとうございました。

• 次のような同様の原則である必要があります。

function v=prst1(t,y)
global m b k Om D F omeg
v(1)=....;
v(2)=y(1);
v=v(:);

• と：

global m b k Om D F omeg
m=1;
b=10;
k=1000;
F=10;
Om=sqrt(k/m);
omeg=1*Om;
D=b/(2*Om*m);
x0=[0;0];
[t,x]=ode23('prst1',0:0.01:10,x0);
plot(t,x)

• しかし、f(t) と k(t) を取得する方法がわかりません。