大二数学物理方法之理想弦振动--matlab研究二阶偏微分方程

本人大二。数学物理方法课程作业分享，有错误请斧正，谢谢！

下面这张图在振幅上有问题，只作模式参考

 

 

这个是用matlab制作的弦振动演示动图，初速度随横轴的刻度成线性增大 

这个动图是初速度为零的情况 

理想弦振动的规律遵循波动方程，一种物理上常见的二阶偏微分方程。通常采用傅里叶分析将初始条件通过傅里叶级数以无数阶具有正弦函数形式的本征模式的线性叠加进行处理。在实际应用中取有限阶逼近。推荐阅读费曼物理学讲义和3b1b了解有关内容。以下为有关的摘录

部分操作代码在此

clc
clear
n=100;%计算级数1到100
v=10%要展示的级数
c=zeros(100)%描述初位移的一百级参数10*10
L=10;%弦长取值越大取点效果越好，但运算时间差很多
a=0.5;%波速
t=40%时间40刚好为一周期T=2L/a,即两倍弦长比上波速，原理可见费曼物理学讲义,
%可理解为一个小震动在两固定端不断被反射，一个来回就是一个周期
f=cell(v,1)%一行一级数据，共v行
第一种谐波分量算法
for n=1:100
fun1 = @(x) 3/10*L*x.*sin(n*pi*x/L)%由于给定初位移只允许sin模式的谐波存在，故不考虑cos项
fun2 = @(x) 3/20*L*(L-x).*sin(n*pi*x/L)%右半边乘以sin(n*pi*x/L)
fun3 = @(x) x.*sin(n*pi*x/L)%这个是初始速度的函数乘以sin(n*pi*x/L)
c(n) = 2/L*integral(fun1,0,L/3)+2/L*integral(fun2,L/3,L)%函数法求积分得参数，点乘！！！！
d(n) =-(1/pi)*integral(fun3,0,2*L/(a*n))
end
% 第二种谐波分量算法
% L=10;c=zeros(10)
% for n=1:100%数值法求每级系数，精度不高，时间长，吃算力，不推荐
%     for x=linspace(0,1/3*L,500)
%     c(n)=c(n)+3/10*L*x*sin(n*pi*x/L)*0.01
%     end
%     for x=linspace(1/3*L,L,1000)%加上0.01是为了避免重复
%         c(n)=c(n)+3/20*L*(L-x)*sin(n*pi*x/L)*0.01
%     end
% end
% 第三种谐波分量算法：复数形式的傅里叶系数取实数复数，由于本初始条件特殊，不列。


%%
hold on
for n=1:v
x=0:0.01:L
%figure(n)
f{n,1}=(c(n)*sin(n*pi/L*x))*cos(t*a*n*pi/L)+(d(n)*cos(t*a*n*pi/L))*sin(n*pi/L*x);%求出每一行所在级的一列初位移
subplot(2,5,n)
plot(x,f{n,1})%画出每一级的初位移分量
title(['Stage ',num2str(n)])
xlabel('length')
end
hold off
sum(numel(x))=0;
for i=1:numel(x)
    summ=0
    for j=1:v
        summ=summ+f{j,1}(i)
    end
    sum(i)=summ
end
figure(11)
plot(x,sum)
title(['time=',num2str(t),' stages added=',num2str(v)])

下面avi视频生成代码是在csdn找到的，感谢分享。与以上代码部分重复

%create files for string oscillation
clc
clear
n=100;%计算级数1到100
c=zeros(10)%描述初位移的一百级参数10*10
L=10;%弦长取值越大取点效果越好，但运算时间差很多
a=0.5;%波速
%t为时间
f=cell(n,1)%一行一级数据，共v行
% for n=1:100
% fun1 = @(x) 3/10*L*x.*sin(n*pi*x/L)%左半边
% fun2 = @(x) 3/20*L*(L-x).*sin(n*pi*x/L)%右半边
% fun3 = @(x) x.*sin(n*pi*x/L)%这个是初始速度的函数乘以sin(n*pi*x/L)
% c(n) = 2/L*integral(fun1,0,L/3)+2/L*integral(fun2,L/3,L)%函数法求积分得参数，点乘！！！！
% d(n) =2/(n*pi*a)*integral(fun3,0,L)
% end
for n=1:100
fun1 = @(x) 3/10*L*x.*sin(n*pi*x/L)%由于给定初位移只允许sin模式的谐波存在，故不考虑cos项
fun2 = @(x) 3/20*L*(L-x).*sin(n*pi*x/L)%右半边乘以sin(n*pi*x/L)
fun3 = @(x) x.*sin(n*pi*x/L)%这个是初始速度的函数乘以sin(n*pi*x/L)
c(n) = 2/L*integral(fun1,0,L/3)+2/L*integral(fun2,L/3,L)%函数法求积分得参数，点乘！！！！
d(n) =-(1/pi)*integral(fun3,0,2*L/(a*n))
end



%counter=0%为了配合avimaker的文件命名规则
% for t=0:40
%     counter=counter+1
%   for n=1:100
%   x=0:0.01:L
%   f{n,1}=c(n)*sin((n*pi)/L*x)*cos(t*a*n*pi/L)%求出每一行所在级的一列初位移
%   end
% sum(numel(x))=0;%sum为100级的叠加总位移
%   for i=1:numel(x)
%     summ=0
%       for j=1:100
%         summ=summ+f{j,1}(i)
%       end
%     sum(i)=summ
%   end
  %save .txt -ascii sum%生成sum向量的txt文件在系统路径 //eval(['save ','stringmov',num2str(counter) '.txt sum -ASCII -double'])%生成文件名
  %eval(['save ','stringmov',num2str(counter) '.txt sum -ASCII '])%生成文件名
%end
t = 0:40;   % 设置间隔
% 设置保存对象
aviobj = VideoWriter('stringmov2.avi');   % 创建一个名为***.avi的AVI文件
aviobj.FrameRate = 8;   % 画图的帧率，越大表示越快10
open(aviobj);   % 打开AVI文件

for t=0:40
    
    for n=1:100
        x=0:0.01:L
        f{n,1}=sin(n*pi/L*x)*(c(n)*cos(t*a*n*pi/L)+d(n)*sin(t*a*n*pi/L))%求出每一行所在级的一列初位移
    end
    sum(numel(x))=0;%sum为100级的叠加总位移
    for i=1:numel(x)
        summ=0
        for j=1:100
            summ=summ+f{j,1}(i)
        end
        sum(i)=summ
    end
    plot(x,sum)
    ylim([-50 50])
    title(['time=',num2str(t)])
    xlabel('x')
    ylabel('u')
    CurrFrame = getframe;   % 获取每帧图片
    writeVideo(aviobj, CurrFrame);   % 写入AVI文件中
end
% 关闭AVI文件
close(aviobj);