数学物理方程之基于数值计算方法的弦振动方程求解2数学物理方法中的平行四边形法则目录摘要、关键词…………………………………………… 2页有限差分法介绍………………………………………… 3页程序描述………………………………………………… 6页计算机处理……………………………………………… 8页Matlab作图…………………………………………… 10页特别鸣谢………………………………………………… 11页摘要、关键词摘要：继上次关于弦振动方程的“平行四边形法则”求解之后，我们又从数值计算的角度入手，对弦振动方程进行计算和模拟，从而验证“平行四边形法则”解弦振动方程的正确性。

关键词：有限差分法、数值计算、弦振动方程附： 弦振动方程：4(0,)(1,)0(,0)(1),(,0)8tt xxtu uu t u tu x x x u t x =⎧⎪==⎨⎪=-=⎩211((1))()'()()''()()+()()2!n n nu i h u ih u ih h u ih h u ih h -=+-+-+-……！211((1))()'()''()+()2!n n nu i h u ih u ih h u ih h u ih h +=+++ ……！()((1))'()()u ih u i h u ih o h h--=+((1))()'()()u i h u ih u ih o h h+-=+2((1))((1))'()()2u i h u i h u ih o h h+--=+有限差分法介绍以弦振动方程为例：2(,)(0,)(,)0(,0)()(,0)()tt xx t u a u f x t u t u l t u x x u x x ⎧=+⎪==⎪⎨=Φ⎪⎪=ψ⎩对于一定的u （x ，t ），我们用“差分”代替“微商”，从而将 数差值描述，可得：以及将第一个式子的右边第一项移至左边,得： ^…同理可得， 两式做差：22((1))((1))ih =h u i h u i h u +--（）(,)(,)ni u x t u i x n t u =∆∆=1122(,)n n n i i i tt tt u u u u u i n t +--+==∆1122(,)n n n i i i xx xx u u u u u i n x +--+==∆21122(,)n n n i i i tt uu u u a f i n x+--+=+∆2222ta r x∆=∆ 2122122112(1)(,)n n n n n i i i i iu r u r u r u n t f i n ---+-=+-+-+∆用中心差分的一阶导数表示二阶导数，化简： 由此引入 则 则弦振动方程 可以表示为：我们定义 为网格比则由此可知，每一个格点u （i ，11(,0)()()2i it u u u x x i t--=ψ=ψ=∆(,0)t u i 1i u 202020221121221221100,/10.5(2(1)2()(,)0,0/12(1)(,)ni i i i n n n n i i i i i l x u r u r u r u t i x t f i x n t n i l x r u r u r u n t f i x n t +----+-=∆-⎧⎪=+-++∆Φ∆+∆∆∆=<<∆-⎨⎪+-+-+∆∆∆⎩ 其他n)均由u （i+1，n+1）、u （i ，n ）、u(i-1,n-1)、 u(i,n-2)等其余四点所确定：由此我们可以采用“递归”的思想，借助计算机进行快速计算，从而得到各个格点的值.值得注意的是,①在边界上u ≡0.②在初始层上的点(即u （i ，0））无法用上述公式计算，还需借助初始条件，即：012020201211(,0)()2(1)(,1)i i i i i i u i u x u r u r u r u u t f i -+-∴==Φ=+-+-+由 和 两式相加，消去可得020*********.5(2(1)2()(,)i i i i u r u r u r u t i x t f i x n t +-=+-++∆Φ∆+∆∆∆综上：届此，我们可以将此式编入程序（采用“递归”思想），详细代码见下一节。

程序描述C语言：#include <stdio.h>#include <stdlib.h>#define imax 10 //定义长、宽格点数，其数值参考计算机性能而定#define jmax 60double f(double i,double n);double u(int i,int n);double FEI(double x);double PUSI(double x);double l=1,t=3; //设l=1，t=3double dert_t,dert_x; //长宽格距double R;int main (){int i,n,a;double A[imax][jmax];FILE *fp;dert_x=l/imax;dert_t=t/jmax;a=2; //a的值R=a*a*dert_t*dert_t/dert_x/dert_x;printf("%lf\n",R);if((fp=fopen("G:\\matlab\\myfile.txt","r+"))==NULL) //创建文件{printf("cannnot open file!\n");return 0;}for(n=0;n<jmax;n++){for (i=0;i<imax;i++){A[i][n]=u(i,n);printf("%9.5lf ",A[i][n]);fprintf(fp,"%9.5lf",A[i][n]);fputs(" ",fp);}printf("\n");fputs("\n",fp);}fclose(fp);return 0;}double u(int i,int n) //函数u(x,t){double result;if(i==0||i==imax-1){result=0;}else{if(n==0){result=FEI(i*dert_x);}else if(n==1){result=0.5*(R*R*u(i+1,0)+2*(1-R*R)*u(i,0)+R*R*u(i-1,0)+2*dert_t*PUSI( x)+dert_t*dert_t*f(i*dert_x,n*dert_t));//result=0.5*(u(i+1,0)+u(i-1,0)+2*dert_t*PUSI(i*dert_x));}else{result=R*R*u(i+1,n-1)+2*(1-R*R)*u(i,n-1)+R*R*u(i-1,n-1)-u(i,n-2)+ dert_t*dert_t*f(i*dert_x,n*dert_t);//result=u(i+1,n-1)+u(i-1,n-1)-u(i,n-2);}}return result;}(2)nO double FEI(double x) //u(x,0)=FEI(x)初始条件 {double fei; fei=x*(1-x); return fei; }double f(double i,double n) //外强迫项f(x,t) {double f;f=0; //此函数可以由题变化，此处暂定为0 return f; }double PUSI(double x) //Ut(x,t)=PUSI(x,t) {double PUSI;PUSI=8*x; //一阶导数 return PUSI; }计算机处理在实际操作过程中，鉴于递归调用在后期计算时间上的复杂度以及个人电脑（PC ）的性能，我们仅对上节课中所涉及到的弦振动方程（第二页所附）进行了模拟，故对程序中的函数关系式进行了修改（即采用了“//”后面的关系式）：注：有程序代码可知，第一行为r2的值：1.000000，随后的数组均输入到文件“myfile.txt”中.文件中数据如下：随后我们可以借助matlab对该文档进行调用，从而绘图，详见下一节。

Matlab作图代码如下（附截图）：M=load('G:\matlab\myfile.txt')x=(0:0.9:0.1)t=(0:2.95:0.05)surf(M)xlabel('x')ylabel('t')zlabel('u(x,t)')title('弦振动方程数值计算方法-surf')图像：或。