二阶微分方程画三维图的问题

lz2239129

设计要求画三维图，顿时崩溃。不会编程，求大神们指教！



其中eta1包含B0z，要求B0z是x轴。画出B0z x 和的三维图。
我编的程序是
方程文件：
function dx=Yqiangpo_F(t,x,flag,eta2,P0,omega,eta1,omegan)
dx=[x(2);-eta1*x(2)-omegan^2*x(1)+eta2*P0*sin(omega*t)];


求解画图：
clear
clc
Omega=2*pi*150;
n=60*Omega/(2*pi);
rho=2670;
sigma0=3.63*(10^7);
mu=0.34;
E=71*(10^9);
h=0.01;
P0=5000;
B0z=0:0.04:4;
R=0.5;
DM=E*h.^3./(12*(1-mu^2));
omega=2*pi*50;
C1=-(6+2*mu)/(5+mu); C2=(1+mu)/(5+mu)
A1=rho.*h.*R.^2/2+((rho.*h.*R.^2)/2-(rho*h.^3)/6)*C1+(rho*h.*R.^2/3-rho*h.^3/3)*C2+(rho*h.*R.^2/6-rho*h.^3/12)*C1^2+(rho*h.*R.^2/10-rho*h.^3/6)/C2^2+(rho.*h.*R.^2/4-5*rho*h.^3/18)*C1*C2;
A2=-sigma0*h.^3.*B0z.^2.*(2*C1+4*C2+C1^2+C2^2+10*C1*C2/3)/12;
A3=DM.*(32*C2^2/3+12*C1*C2+32*C2)./R.^2+rho*h.^3*Omega^2*(2*C1+4*C2+C1^2+2*C2^2+10*C1*C2/3)/12;
A5=(1/2+(C1)/4+(C2)/6).*R.^2;
omegan=sqrt(A3./A1);
qst=(A5*P0)./(A1.*omegan.^2)
s=omega./omegan;
epsilon=A2./(2*omegan.*A1);
eta2=A5/A1;
eta1=A2/A1;

beta=1./sqrt((1-s.^2).^2+(2*epsilon.*s).^2);

x0=[0.00001,0.1];
t=0:0.1:10;
options=odeset('RelTol',1e-6,'AbsTol',1e-6);
[t,x]=ode45('Yqiangpo_F',t,x0,options,eta2,P0,omega,eta1,omegan);
[B,X1]=meshgrid(B0z,x(1));
surf(B,X1,x(2))



程序前面大部分都是系数。我想画B x1 x2 的关系。但是图中ode45开始就出错了。请问大家应该怎么改？

lz2239129

不要沉啊

lz2239129

大神都不在嘛:Q

lz2239129

:lol:lol:lol:lol:lol:lol:lol:lol:lol

lin2009

ode45等函数的是求解一阶微分方程组的。
因此首先要先把高阶的微分方程转化为一阶方程组的形式再求解，可能还需要用到实际电路的其它条件。

另外，若eta1，eta2为常系数的话，这个二阶微分方程应该比较容易求解，即使是手工求解。
很显然，解的形式应为： x = c1* sin(omega *t) + c2*cos(omega*t)

如果eta1，eta2不为常系数，则需转化为一阶微分方程组。