用MATLAB仿真下图的公式，参数什么的我就不写了，直接粘贴我的程序，新手，求教

你复制以下代码到你的电脑里，然后按你的要求修改。

clear all,close all,clc,clf

t=0:1e9:10e9;   %可以改

fn=length(t);

i=1;

while i<=fn

T=2175e-12;

t0=t(i);

H1=cos(145015501550e-13t0t0/12+pi/2);

syms k;

n=100;    %明确n次方值

f=exp(j2t(k-1)T);

H2=symsum(f,k,1,n);

H=H1H2;

fudu0=abs(sum(H));

fudu1=double(fudu0);

fudu(i)=fudu1;

i=i+1;

end

fudu1=10log10(fudu/max(fudu));

plot(t,fudu1,'r-')

1、用数值方法求解微分方程（未指定初始条件，按零初始条件考虑）：

% 常数定义

m=196;

k=19600;

c=2940;

clf

tstr = {'忽略阻尼', '考虑阻尼'};

for n=1:2

    subplot(2,1,n)

    dx=@(t,x)[x(2); (160sin(19t)-kx(1)-(n-1)cx(2))/m];

    [t,x]=ode45(dx,[0 3],[0 0]);

    ax=plotyy(t,x(:,1),t,160sin(19t));

    %legend([h1 h2],'x(t)','P(t)')

    xlabel t

    axes(ax(1)); ylabel x(t)

    axes(ax(2)); ylabel P(t)

    title(tstr{n})

end

对于考虑阻尼影响的情况，系统稳态响应为和输入相同频率的正弦波，由图中的峰值可以大致计算出幅值的放大倍数以及相位滞后。而对于忽略阻尼影响的情况，由于输出由自由振荡和强迫振荡两部分组成，且频率不同，呈现出来的响应曲线不太容易辨别幅值放大及相位变化。事实上，无阻尼系统在传统控制理论中认为是不稳定的，那么建立在稳态响应基础上的频域方法从概念上说也是有疑问的。

2、理论方法求解。

考虑两种做法，一是求解微分方程的解析解：

>> x1=dsolve('196D2x+19600x=160sin(19t)','x(0)=0,Dx(0)=0')

x1 =

76/12789sin(10t)-40/12789sin(19t)

>> x2=dsolve('196D2x+2940Dx+19600x=160sin(19t)','x(0)=0,Dx(0)=0')

x2 =

6308/2845871exp(-15/2t)sin(5/27^(1/2)t)7^(1/2)+1900/1219659exp(-15/2t)cos(5/27^(1/2)t)-580/406553sin(19t)-1900/1219659cos(19t)

其中x1为无阻尼情况，包括两项，其中后面一项为强迫振荡，可知放大倍数为40/12789/160 = 19548e-005，相角滞后为0；x2为考虑阻尼的情况，结果包括4项，前两项为衰减项，稳态响应只有后两项，则其幅值放大倍数和相角滞后分别为

>> norm([-1900/1219659, -580/406553])/160

ans =

  13202e-005

>> atan2(-1900/1219659, -580/406553)/pi180

ans =

 -1324831

二是用频率特性的概念直接求：

G = tf(1,[m c k]);

w = 19;

Gjw = evalfr(G,jw);

mag = abs(Gjw)

phase = angle(Gjw)180/pi

得到的

mag =

  13202e-005

phase =

 -1324831

即分别为幅值放大倍数和相角滞后，和上面求解微分方程的结果一致。

问题

首先，说一下条件给的不清楚或可能有错的地方：

1、按照公式中使用的符号，下面的常数d、h应为D、H，a1、a2应为R1、R2；

2、公式中，第三重积分的下限确定是2sqrt(R1^2-x)？后面有讨论

3、求积分得到的结果就是一个具体的数，标题要求的绘制仿真图是什么意思？

参考代码

下面给出我写的代码，供参考。

% 常数定义 

I=1; u0=4pi1e-007; D=015; H=114; R1=75; R2=14; 

  

% 六元函数定义（注意把里面第二重积分的系数以及D^3折算到函数中） 

f=@(t,h,r,y,x,z)(-rsin(t)(y-sin(t))-rcos(t)(x-cos(t))) /  

    ((x-rcos(t))^2+(y-rsin(t))^2+(z-h)^2)^(3/2)  u0I^2/(4pi)/D^3; 

  

% 积分参数设置 

tol = 1e-6; 

quadf = @quadl;

trace = 1;

  

% 先使用triplequad计算内三重积分（积分限均为常数，h、r注意换算） 

Ithr = @(y,x,z) triplequad(f,0,2pi,0,DH,DR1,DR2,tol,quadf,y,x,z); 

  

% 由于被积分函数必须能够接受向量输入并返回向量输出，所以外面的三重积分必须逐层进行

Iy = @(x,z) quadf( @(Y)arrayfun(@(y)Ithr(y,x,z),Y),2sqrt(R1^2-x),sqrt(R2^2-x^2),tol,trace);

Ix = @(z) quadf( @(X)arrayfun(@(x)Iy(x,z),X),R1,R2,tol,trace);

Iz = quadf( @(Z)arrayfun(@(z)Ix(z),Z),0,DH,tol,trace);

注解及讨论

1、对代码自身不想再多做解释，该说的注释里都说了。自己好好看看，实在看不懂再问吧（要更好地理解代码，需要好好熟悉一下匿名函数）。

2、需要注意：重积分的计算时间不是根据积分的重数成倍增加，而是按照幂级数规律增加的，所以，计算六重积分所需的时间是比较长的（这段代码在我电脑上求解时间近8分钟）。

3、按照当前的程序，求出的结果：

Iz =

  11554e-005

由于y的积分下限看起来有点怪，像是sqrt(R1^2-x^2)之误，但仔细分析应该没错。如果积分下限取sqrt(R1^2-x^2)，则计算结果为复数：

Iz =

  54730e-005 -67361e-005i

原因是，对x的积分范围是R1-R2，也就意味着，x的取值范围最大可以是R2，而里面那一重积分的下限是sqrt(R1^2-x^2)，当x取值大于R1时就会出现复数。

4、上面的代码侧重于实现方法，具体的数学表达式不保证绝对没问题，请自行仔细核实。

如果x未知或者你需要把它作为一个变量，用符号运算，syms。分别赋值被加函数至一个符号数组，然后求和。

或者给x赋具体数值（一般是一定范围内的一个数列）后，同理赋值被加函数至一个数组，求和。

>

以上就是关于用MATLAB仿真下图的公式，参数什么的我就不写了，直接粘贴我的程序，新手，求教全部的内容，包括:用MATLAB仿真下图的公式，参数什么的我就不写了，直接粘贴我的程序，新手，求教、怎么用matlab编有阻尼d簧系统的仿真程序，题3.1、求用matlab计算一个六重积分的具体程序，绘制仿真图等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！