matlab二阶节的系数,如何利用matlab求解矩阵系数的二阶微分方程
M=[2,0;0 1 ]; %质量矩阵
K=[6 -2;-2 4]; %刚度矩阵
a=0;b=0;
C=a*K+b*M;
dt=0.28;
t=0:dt:2.8;
ft0=zeros(length(K),length(t));
for i=1:length(t)
ft0(1,i)=10; %在节点4的竖直方向加大小为200N的阶跃力
end
dsp=zeros(length(K),length(t)); % 位移
vel=zeros(length(K),length(t)); % 速度
acc=zeros(length(K),length(t)); % 加速度
%--------------------------------------------------------------------------
% (2) Newmark
alpha=0.3; beta=0.6; % 稳定条件
acc(:,1)=inv(M)*(ft0(:,1)-K*dsp(:,1)-C*vel(:,1)); % 计算初始加速度 (t=0)
ekk=K+M/(alpha*dt^2)+C*beta/(alpha*dt); % 计算有效刚度矩阵
for it=2:length(t)-1 % 时间步循环
cfm=dsp(:,it)/(alpha*dt^2)+vel(:,it)/(alpha*dt)+acc(:,it)*(0.5/alpha-1);
cfc=dsp(:,it)*beta/(alpha*dt)+vel(:,it)*(beta/alpha-1)...
+acc(:,it)*(0.5*beta/alpha-1)*dt;
efd=ft0(:,it)+M*cfm+C*cfc; % 计算有效力矢量
dsp(:,it+1)=inv(ekk)*efd; % t+dt时刻的位移
acc(:,it+1)=(dsp(:,it+1)-dsp(:,it))/(alpha*dt^2)-vel(:,it)/(alpha*dt)...
-acc(:,it)*(0.5/alpha-1); % t+dt时刻的加速度
vel(:,it+1)=vel(:,it)+acc(:,it)*(1-beta)*dt+acc(:,it+1)*beta*dt; % t+dt时刻的速度
end
plot(t, dsp(2,:),'-b*')
hold on
xlabel('Time(seconds)')
ylabel(' Vertical displ. (m)')
% (2) Wilson
%--------------------------------------------------------------------------
theta=1.4; % 稳定条件参数
acc(:,1)=inv(M)*(ft0(:,1)-K*dsp(:,1)-C*vel(:,1)); % 计算初始加速度 (t=0)
ekk=K+M*(6/(theta*dt)^2)+C*(3/(theta*dt)); % 计算有效刚度矩阵
for it=2:length(t)-1
cfm=dsp(:,it)*(6/(theta*dt)^2)+vel(:,it)*(6/(theta*dt))+2*acc(:,it);
cfc=dsp(:,it)*(3/(theta*dt))+2*vel(:,it)+acc(:,it)*(theta*dt/2);
efd=ft0(:,it)+theta*(ft0(:,it+1)-ft0(:,it))+M*cfm+C*cfc; %计算有效力矢量
dtheta=inv(ekk)*efd; % t+ dt时刻的位移
acc(:,it+1)=(dtheta-dsp(:,it))*(6/(theta^3*dt^2))...
-vel(:,it)*(6/(theta^2*dt))+acc(:,it)*(1-3/theta); % t+dt时刻的加速度
vel(:,it+1)=vel(:,it)+acc(:,it+1)*dt/2+acc(:,it)*dt/2; % t+dt时刻的速度
dsp(:,it+1)=dsp(:,it)+vel(:,it)*dt...
+(acc(:,it+1)+2*acc(:,it))*(dt^2/6); % t+dt时刻的位移
end
plot(t, dsp(2,:),'-g*')
hold on
xlabel('Time(seconds)')
ylabel(' Vertical displ. (m)')
%--------------------------------------------------------------------------
% (2) 振型叠加法
%施加不平衡激振力
t=t';
[V,D]=eig(K,M); % 计算特征值和特征向量
[lambda,ki]=sort(diag(D)); % 给特征值和特征向量排序
omega=sqrt(lambda); % 角频率w.
omega1=sqrt(lambda)/(2*pi); % 固有频率 Hz.
V1=V(:,ki);
h=V1'*M*V1;
Factor=diag(V1'*M*V1); %模态质量
Vnorm=V1*inv(sqrt(diag(Factor))); % 正则化振型向量
VnormK=diag(Vnorm'*K*Vnorm); %模态刚度
omega0=diag(sqrt(Vnorm'*K*Vnorm)) ; % 正则化模态刚度
VnormM=diag(Vnorm'*M*Vnorm); %正则化模态质量
Fnorm=Vnorm'*ft0 ; % 模态力矢量
Modamp=Vnorm'*(a*M+b*K)*Vnorm; % 模态阻尼矩阵
zeta=diag((1/2)*Modamp*inv(diag(omega0))); % 阻尼比
%--------------------------------------------------------------------------
% % (4) 模态坐标的响应
%--------------------------------------------------------------------------
q0=zeros(length(K),1);
dq0=zeros(length(K),1); % 初始化位移和速度
eta0=Vnorm'*M*q0; deta0=Vnorm'*M*dq0; %2.2-53% 初始条件模态坐标的位移和速度
eta=zeros(length(t),length(K)); % %初始化位移2.2-52
for i=1:length(K) % t(i)时刻的响应
phase0=omega0(i)*t; %w*t omega0为无阻尼固有频率
omegad=omega(i)*sqrt(1-zeta(i)^2);
phase=omegad*t; %wd*t其中(omegad为有阻尼固有频率)
Exx=exp(-zeta(i)*omega(i)*t); %中间变量
C1=eta0(i); %初始位移
C2=(deta0(i)+eta0(i)*zeta(i)*omega0(i))/omegad; %中间变量
D1=zeta(i)*omega(i)/omegad;
II=ones(length(t),1);
XX=Fnorm(i)/(omegad^2+zeta(i)^2*omega(i)^2);
eta(:,i)=C1*Exx.*cos(phase)+C2*Exx.*sin(phase)+ XX*(II-Exx.*cos(phase)-D1*Exx.*sin(phase));%有阻尼系统
% eta(:,i)=Fnorm(i)/((omega(i))^2)*(1-cos(omega(i)*t)); %无阻尼系统
end
%--------------------------------------------------------------------------
% (6) 将模态坐标转化到物理坐标系
%--------------------------------------------------------------------------
eta=eta'; %模态坐标
y=Vnorm*eta; %物理坐标
plot(t, y(2,:),'-r*')
xlabel('Time(seconds)')
ylabel(' Vertical displ. (m)')
matlab二阶节的系数,如何利用matlab求解矩阵系数的二阶微分方程相关推荐
- 利用MATLAB进行系统时域分析,实验二 利用matlab进行系统的时域分析
实验二 利用matlab进行系统的时域分析 实验二 利用MATLAB进行系统的时域分析 1.实验目的 在理论学习的基础上,通过本实验熟悉LTI连续时间系统的时域分析方法, 熟悉系统的零输入响应.零状态 ...
- matlab系列之(一)——利用matlab实现任意两个多项式相加
上课时的课程作业,后续我会持续整理出来注释好,供大家共同学习!致谢课程老师! 一.问题描述 输入任意两个多项式,相加后输出结果: 二.问题分析 输入多项式可以采用输入系数矩阵或完整多项式的方式,为了符 ...
- 用matlab编程实现h鲁棒控制算法,利用matlab实现H-infinity鲁棒控制
利用matlab实现H-infinity鲁棒控制 利用Matlab实现H∞控制 Prof. Dr.-Ing.F.Allgwer Institute for Systems Theory and Aut ...
- 在matlab中实现累乘,如何利用matlab设计一个线性相位FIR带通滤波器,并在FPGA上实现...
设计要求 利用matlab设计一个线性相位FIR带通滤波器,并在FPGA上实现. 1.滤波器指标:过渡带带宽分别为100~300HZ,500~700HZ,阻带允许误差为0.02,通带允许误差为0.01 ...
- matlab求多元函数的极小值,[转载]利用MATLAB求多元函数的极值(2)
利用MATLAB求多元函数的极值分两种情况,(1)无约束条件:(2)有约束条件. (2)有约束条件下求极小值的方法: 假设多变量非线性函数的数学模型为 min f(x) c(x)<=0 ceq( ...
- matlab抓取网页信息,如何利用Matlab抓取网页数据
如何利用Matlab抓取网页数据 2019-01-01 %朋友需要做金融方面的分析,要求从网站上下载大量的数据,一个一个复制粘贴太费事.我写了一个简单的网络爬虫,主要用到正则表达式,可以自动下载网页源 ...
- matlab 随机相位的正弦信号,利用MATLAB绘制随机相位正弦波.docx
实验二 利用MATLAB绘制随机相位正弦波的均值,方差和自相关函数的图像[实验目的]通过绘制图像,深入理解随机相位正弦波的均值,方差和自相关函数.[实验学时]课外完成[实验准备]1.熟悉随机相位正弦波 ...
- matlab 画三条曲线,如何利用MATLAB(plot 3函数和fplot3函数)绘制三维曲线?
文章目录 0 前言 1 plot3函数 1.1 plot3函数的基本用法 1.2 plot3(x,y,z)函数参数的变化形式 1.3 含多组输入参数的plot3函数 1.4 含选项的plot3函数 2 ...
- 利用matlab求零输入响应波形,实验3 利用matlab求LTI连续系统的响应
实验3 利用matlab求LTI连续系统的响应 一. 实验目的: 1. 了解LTI系统的冲激响应h(t)及matlab实现: 2. 了解LTI系统的阶跃响应g(t)及matlab实现: 3. 了解LT ...
- matlab元胞自动机学风演化,利用MATLAB和VC60混合编程技术研究元胞自动机动态演化过程...
利用MATLAB和VC60混合编程技术研究元胞自动机动态演化过程 第! !卷!第期 ! 成都理工大学学报! 自然科学版 !# $ % - 2: 5 34- 1 -6;!9 ? : A ? $ % $ ...
最新文章
- 企业内网中的WSUS更新服务 服务器连接到Microsoft Update来获取更新程序
- C++ Primer 5th笔记(chap 14 重载运算和类型转换)可调用对象与function
- 微型计算机2018读者调查,《微型计算机》2018年度电竞品牌影响力调查获奖读者揭晓!...
- Codeforces Round #564 (Div. 2)A
- P3538-[POI2012]OKR-A Horrible Poem【hash,字符串】
- 计算机网络学习笔记(9. 报文交换与分组交换③)
- AMP (LAMP/WAMP)
- android+JPEG+编码,Android_解析:android 如何从JPEG生成BufferedImage,如下所示:复制代码 代码如下 - phpStudy...
- 性能测试中容易混淆的概念
- 移动目录下的隐藏文件
- java程序设计基础_陈国君版第五版_第十章习题
- OA协同办公软件评测 —— Tower篇
- MineCraft建模工具
- win10 u盘 修复计算机,U盘启动盘修复win10系统的方法
- 电脑计算机管理摄像头服务,电脑上打开摄像头的方法
- html 播放微信amr音频文件,如何在微信中播放amr格式的文件?
- TQ2440内核linux2.6.28移植
- 基频和倍频的概念_基频峰,泛频峰,倍频峰,二倍频峰的区别
- win10家庭版不能远程连接,升级企业版过程
- 从0开始的编程学习计划
热门文章
- Windows下UDP编程
- hsqldb mysql 语法_[spring batch]建表语句(hsqldb改mysql)
- 郑州大学c语言实验报告答案,郑州大学c语言实验报告册答案
- 专利 | word图片设置为黑白
- python mql4跟单_MT4多功能本地跟单EA
- 大乐透号码随机生成(仅供参考学习)
- NRF52840学习历程(四)定时器
- 计算机网络应用模拟试卷,《计算机网络应用基础》模拟试卷(八)(附答案)
- DataGear 国产开源BI
- 数据结构-串、数组和广义表