MATLAB解微分方程

来源:百度知道 编辑:UC知道 时间:2024/05/21 23:40:43
4x”+35.78x’+2000x=F(t)
F(t)=100 0<t<=0.1
F(t)=-1000(t-0.2) 0.1<t<=0.2
要求用ode45解出x'和x(数值解)
F(t)用如下的function表示
function [force]=f(t)
if t<0.1
force=100
elseif t<0.2
force=-1000*(t-0.2)
else
force=0
end
谢谢了
t=0时刻,x=0;x'=0

终于解出来了,解析解与数值解不同,解析解正确,数值解问题出在t=0.1时的初值选取上。

function hh
clc;clear
x0=[0 0];
[T1,X1] = ode45(@odefun1,0:0.01:0.1,x0)
[T2,X2] = ode45(@odefun2,0.1:0.01:0.2,[X1(end,1),X1(end,2)])
plot(T1,X1(:,2),T2,X2(:,2))
grid
%解析解
x1=dsolve('4*D2x+35.78*Dx+2000*x=100','x(0)=0','Dx(0)=0')
x2=dsolve('4*D2x+35.78*Dx+2000*x=-1000*(t-0.2)','x(0)=0','Dx(0)=0')
t1=0:0.01:0.1;
t2=0.1:0.01:0.2;
x1=subs(x1,'t',t1)
x2=subs(x2,'t',t2)
figure
plot(t1,x1,t2,x2)
grid
function dx=odefun1(t,x)
x1=x(1);x2=x(2);
dx1=x1;
dx2=(100-35.78*x1-2000*x2)/4;
dx=[dx2;dx1];
function dx=odefun2(t,x)
x1=x(1);x2=x(2);
dx1=x1;
dx2=(-1000*(t-0.2)-35.78*x1-2000*x2)/4;
dx=[dx2;dx1];

结果:
T1 =

0
0.0100
0.0200
0.0300
0.0400
0.0500
0.