Question

MATLAB 求解方程

条件为：t=0时，q0=0,Dq0=0，a=200-2000。
求解q0关于t和a的三维图，然后求解q1-9q0关于时间t和a的三维图，以及a为定值时（a=1500）的q1-9q0的曲线图。

Answer 1

以前回答过几个相关的问题，大多数代码的形式都是用过的，不多作解释了。

需要说明一下的是，最后一项的导数需要首先展开成 2*(q1-9*q0)*(dq1-9*dq0)。

参考代码：

syms t
q1=0.005873-0.0002634*cos(42.41*t)-0.0003364*cos(84.82*t);
dq1=eval(['@ (t)' char(diff(q1))]);
q1=eval(['@ (t)' char(q1)]);
T=linspace(0,0.5,1000);
N = 5;
A=linspace(200,2000,N);
X=zeros(1000,N);
for k=1:N
    a=A(k);
    dx=@ (t,x) [x(2); ...
        1444408.2*2*(q1(t)-9*x(1))*(dq1(t)-9*x(2))-535.81*a*x(1)-0.5*x(2)];
    [t,x]=ode45(dx,T,[0 0]);
    X(:,k)=x(:,1);
end

mesh(A,T,X)
xlabel a; ylabel t; zlabel q_0

上面的代码是画出q0随a及t变化的曲面（为节省计算时间，t的范围取0-0.5，a只取了5个值）：

q1随时间变化的曲线可以直接画出：

plot(T,q1(t))

如果把q1和9*q0画在一起，可见二者变化幅度基本相当，但有一个较大的移位：

plot(T,q1(t),T,X*9)
xlabel t; ylabel 'q_1, 9*q_0'

如果把q1的常数项减去：

plot(T,q1(t)-0.005873,T,X*9)
xlabel t; ylabel 'q_1-0.005873, 9*q_0'

按照题主要求的q1-9*q0画曲面：

mesh(A,T,repmat(q1(T(:)),1,N)-9*X)
xlabel a; ylabel t; zlabel q_1-q_0

至于说a为定值时（a=1500）的曲线，给a赋值，然后很容易做，不赘述了。