微分方程数值解课程设计题目
微分方程数值解法C语言-课程设计
微分方程数值解法C语言-课程设计微分方程数值解法C语言由于对matlab语言不熟悉，所以还是采用C。
前面几个都比较简单，最后一个需要解非其次方程组。
采用高斯—Jordan消元法（数值分析）求逆解方程组，也再一次体会到算法本身的重要性，而不是语言。
当然，矩阵求逆的算法也在100个经典的C语言算法之列。
不过偏微分方程数值解的内容的确比较高深，我只能停留在编这种低级的东西的自娱自乐中。
不过解决计算机、数学、信计专业的课程设计还是足够了。
由于篇幅所限，只把源代码粘贴在这。
一：预报矫正格式#include <math.h>#include<iostream>#include<stdlib.h>double count_0( double xn,double yn){//矫正格式double s;s=yn+0.1*(yn/xn*0.5+xn*xn/yn*0.5);return s;}double count_1(double xn,double yn,double y0){//预报格式double s;s=yn+0.05*((yn/xn*0.5+xn*xn/yn*0.5)+(y0/xn*0.5+xn*xn/y0*0.5));return s;}void main(){//计算,步长为0.1,进行10次计算,设初始值double xn=1,yn=1;int i=1;while(i<=10){printf("%16f ,%1.16f ,%1.16f\n",xn,yn,count_1(xn,yn,count_0(xn,yn)));xn=xn+0.1;yn=count_1(xn,yn,count_0(xn,yn));i++;}}二显示差分格式#include<iostream>#include<math.h>#include<stdlib.h>main(){double a[6][11];//初始化；for(int i=0;i<=5;i++){a[0]=0;a[10]=0;}for(int j=1;j<10;j++){double p=3.14*j*0.1;a[0][j]=sin(p);}//按显示格式计算for(i=1;i<=5;i++)for(j=1;j<10;j++)a[j]=a[i-1][j-1]+a[i-1][j+1]; //输出计算好的矩阵for(i=0;i<=5;i++){for(j=0;j<11;j++)printf("%1.10f ",a[j]);printf("\n");}}三龙阁库塔格式#include <math.h>#include<iostream>#include<stdlib.h>double count_k( double xn,double yn){ double s;s=yn/xn*0.5+xn*xn/yn*0.5;return s;}void main(){//步长为0.1double xn=1,yn=1;int i=1;while(i<=11){printf("%f ,%f\n",xn,yn);double k1=count_k(xn,yn);double k2=count_k(xn+0.05,yn+0.05*k1); double k3=count_k(xn+0.05,yn+0.05*k2); double k4=count_k(xn+0.01,yn+0.1*k3); yn=yn+0.1/6*(k1+2*k2+2*k3+k4);xn=xn+0.1;i++;}}四 CRANK--NICOLSON隐式格式#include<iostream>#include<math.h>#include<stdlib.h>double Surplus(double A[],int m,int n);double * MatrixInver(double A[],int m,int n);double * MatrixOpp(double A[],int m,int n) /*矩阵求逆*/ {int i,j,x,y,k;double *SP=NULL,*AB=NULL,*B=NULL,X,*C;SP=(double *)malloc(m*n*sizeof(double));AB=(double *)malloc(m*n*sizeof(double));B=(double *)malloc(m*n*sizeof(double));X=Surplus(A,m,n);X=1/X;for(i=0;i<m;i++)for(j=0;j<n;j++){for(k=0;k<m*n;k++)B[k]=A[k];{for(x=0;x<n;x++)B[i*n+x]=0;for(y=0;y<m;y++)B[m*y+j]=0;B[i*n+j]=1;SP[i*n+j]=Surplus(B,m,n);AB[i*n+j]=X*SP[i*n+j];}}C=MatrixInver(AB,m,n);return C;}double * MatrixInver(double A[],int m,int n) /*矩阵转置*/ {int i,j;double *B=NULL;B=(double *)malloc(m*n*sizeof(double));for(i=0;i<n;i++)for(j=0;j<m;j++)B[i*m+j]=A[j*n+i];return B;}double Surplus(double A[],int m,int n) /*求矩阵行列式*/ {int i,j,k,p,r;double X,temp=1,temp1=1,s=0,s1=0;if(n==2){for(i=0;i<m;i++)for(j=0;j<n;j++)if((i+j)%2) temp1*=A[i*n+j]; else temp*=A[i*n+j];X=temp-temp1;}else{for(k=0;k<n;k++){for(i=0,j=k;i<m,j<n;i++,j++) temp*=A[i*n+j];if(m-i){for(p=m-i,r=m-1;p>0;p--,r--) temp*=A[r*n+p-1];}s+=temp;temp=1;}for(k=n-1;k>=0;k--){for(i=0,j=k;i<m,j>=0;i++,j--) temp1*=A[i*n+j];if(m-i){for(p=m-1,r=i;r<m;p--,r++) temp1*=A[r*n+p];}s1+=temp1;temp1=1;}X=s-s1;}return X;}void initmat_A(double a[][9],double r){ for(int i=0;i<9;i++)for(int j=0;j<9;j++)a[j]=0;for(i=0;i<9;i++){a=1+r;if(i!=8) a[i+1]=-0.5*r;if(i!=0) a[i-1]=-0.5*r;}}void initmat_B(double b[][9],double r){ for(int i=0;i<9;i++)for(int j=0;j<9;j++)b[j]=0;for( i=0;i<9;i++){b=1-r;if(i!=8) b[i+1]=0.5*r;if(i!=0) b[i-1]=0.5*r;}}void initmat_C(double C[][9]){ for(int i=0;i<9;i++)for(int j=0;j<9;j++)C[j]=0;}void main(){double a[100][11];for(int i=0;i<100;i++)for(int j=0;j<11;j++)a[j]=0;//初始化；for(i=0;i<100;i++){a[0]=0;a[10]=0;}for(int j=1;j<10;j++){double p=4*3.14*j*0.1;a[0][j]=sin(p);}//取h=0.1*3.14,r=0.0005，t=0.0001*3.14*3.14; //得到矩阵a和矩阵bdouble A[9][9];initmat_A(A,0.005);double B[9][9];initmat_B(B,0.005);//B矩阵与Un相乘,en是0;double C[9][9];initmat_C(C);double *A_;A_=MatrixOpp(A[0],9,9);//A矩阵求逆；//A逆*Bfor(i=0;i<9;i++)for(j=0;j<9;j++)for(int s=0;s<9;s++)C[j]+=A_[i*9+s]*B[s][j];//填写a表格for(i=0;i<100;i++){for(j=1;j<10;j++)for(int s=0;s<9;s++)a[i+1][j]+=a[s+1]*C[j-1][s];}//输出表格for(i=0;i<100;i++){for(j=0;j<11;j++)printf("%1.8f ",a[j]);printf("\n");}printf("\n"); printf("\n");//利用精确解，求出表格for(i=0;i<100;i++){for(j=0;j<11;j++)printf("%1.8f",exp(-16*0.0001*0.0005*3.14*3.14*i)*sin(4*j*0.1*3.14));printf("\n");}}。
04实验4微分方程数值解1(4)
0.8000 1.6498 1.6125
return o
0.9000 1.7178 1.6733 x1 x2 x3 x4 1.0000 1.78x48 1.7321
1.4 改进的欧拉公式
例 用改进的欧拉公式解方程
先由向前欧拉公式算出yn+1的预测值 y yny1，2yx再,把y(它0)代 1替梯形
y(1)=y0;
yn1
yn
h[ 2
f
(xn ,
yn )
f
( xn1 ,
re0turn 1.0000 yn001)..]12,00n00000,1, 121,..01985419
for n=1:m
称为改进的欧拉公式，它可写作
k1=feval(f,x(n),y(n));
0.3000 1.2662 0.4000 1.3434
yk 1 yk 1 , k 2h
1,2,, n 1
f
( x0 )
3y0
4 y1 2h
y2
,
f
( xn )
yn2
4 yn1 2h
3yn
后两个公式是由二次插值得到，目的是保持在端点处的精度O(h2) 高阶导数的近似公式一般要用插值多项式得到，下面是二阶
导f数(x的) 近((x似fx0(公axx11式)))((：xx0f(xxa22)) yh0)((x2x1
0.1000 3.9854 0.2924 15.0000 2.0000
0.1500 5.9445 0.6906 15.0000 3.0000
0.2000 7.8515 1.2899 15.0000 4.0000
0.2500 9.6705 2.1178 15.0000 5.0000
《微分方程数值解》课程论文题目1
《微分方程数值解》课程论文题目1、给出Dirichlet 边值问题：()()()(),u x q x u x f x a x b''-+=<< (),()u a u b αβ==和两点边值问题(())()()(),,d d u d u L u p x r x q x u f x a x b d x d x d x =-++=<< (),()u a u b αβ==在均匀网格下的差分格式，编写程序，给出误差阶及误差图。
2、考虑守恒型微分方程（I ）：(())()(),d du Lu p x q x u f x a x b dx dx =-+=<< (),()u a u b αβ==给出其直接差分格式和积分插值法差分格式，采用积分插值法时，数值积分分别利用中矩形公式和梯形公式，编写程序，给出误差阶及误差图。
3、考虑守恒型微分方程（II ）：(())()(),d du Lu p x q x u f x a x b dx dx =-+=<< 01()(),u a u a αα'=+01()().u b u b ββ'=+给出其直接差分格式和积分插值法差分格式，对于边界条件采取直接微分法和积分插值法，编写程序，给出误差阶及误差图。
4、考虑Poisson 方程：(,),(,)u f x y x y G -∆=∈ |(,),u x y αΓ=其中G 是xy 平面上一有界区域，其边界Γ为光滑曲线。
给出方程的五点差分格式和九点差分格式，给出其截断误差，编写程序，给出误差阶及误差图。
5、考虑Poisson 方程：(,),(,)u f x y x y G -∆=∈ |(,),u x y αΓ=其中G 是xy 平面上一圆域、环形域或扇形域，其边界Γ为光滑曲线。
给出方程极坐标形式的差分格式，编写程序，给出误差阶及误差图。
6、考虑Laplace 方程：0,(,),u x y G -∆=∈ |(,),u x y αΓ=其中G 是xy 平面上一有界区域，其边界Γ为光滑曲线。
微分方程数值解课程设(2019年)
微分方程数值解课程设计（五点差分格式测试及溢油事故的扩散行为模拟预测）2019年11月内容：1. 数值测试：椭圆型方程五点差分格式理论结论验证（大作业）包括：（1）方法理论结果验证：断误差分析（讨论是否收敛及收敛阶，通过选取不同的步长，观测步长与误差的关系）。
（2）简单的可视化等。
2. 应用实例: 溢油事故的扩散行为模拟预测船舶溢油事故给海洋环境带来了严重的危害。
溢油事故的发生在危害海洋生态环境的同时，也给当地的渔业及旅游业等造成了不可估量的经济损失。
随着时间的移动，在外界条件的作用下，进入海洋中的溢油可能会发生风化、扩散及漂移等一系列的变化，致使处理它的难度增加。
需要更加高效的溢油回收、处置和清除的技术。
处理溢油的常见方法有微生物分解法、化学处理法以及物理回收法微生物。
分解法主要是依赖于能够对石油进行分解的微生物；化学处理法主要是通过消油剂及特定条件下的燃烧来完成的；物理回收方法主要采用机械装置，如专业收油器等强行回收海面溢油。
从环境保护和资源回收再利用的方面考虑，由于消油剂具有毒性及微生物的分解缓慢，采用机械设备回收海面溢油是最理想的处理方式。
下面我们讨论用围油栏结合收油器装置的处理窄河道溢油问题。
当紧急溢油事故发生时，我们通常采用围油栏和收油装置结合的方法对溢油污染进行紧急处理。
溢油发生时我们先用围油栏对溢油进行拦截使溢油范围缩小到一定程度，再利用收油装置对汇聚起来的溢油进行收集。
然而我们在使用收油器时，需要注意收油器的额定收油量。
在收油的过程中，单位时间收油量应与其额定收油量相同，当比额定收油量大时， 将会导致收油装置的堵塞、失效；当比额定收油量小时，收油装置就不能达到最佳的工作状态。
通过数值模拟，我们可以进一步确定围油栏内的油浓度变化，进而根据收油装置在单位时间内额定收油量来选择合适的收油速度。
当溢油发生后，围油栏会立刻围住溢油，假设河道较窄，因此可认为流动的溢油在Y 方向上的浓度分布是一致的，故只需考虑X 方向上的溢油变化，即该物理问题可通过一维浓度扩散方程来进行模拟（我们采取无量纲的形式，只考虑理想状态下的溢油扩散问题进行模拟）。
微分方程数值解课程设计
课程设计说明书题目：常微分方程数值解法课程设计学院〔系〕：理学院年级专业：计算科学11-1学生：指导教师：教师职称：教授燕山大学课程设计〔论文〕任务书院〔系〕：理学院教学单位：年月日燕山大学课程设计评审意见表摘要本文对常微分方程初值问题现有的数值解法进展了综述研究。
主要讨论了几种常用的数值解法：即欧拉法，改良欧拉法，龙格库塔方法，阿达姆斯PECE,PMECME格式等。
文章最后结合常见数值解法，对较为典型的微分方程模型进展数值求解，探讨了上述数值算法在实际建模问题中的应用。
论文阐述的是常微分方程数值解法的几个问题，通过对以下问题的求解一.比拟Adams四阶PECE模式和PMECME模式。
二.求解贝塞尔方程并与准确解比拟。
三.小型火箭初始重量为1400kg，其中包括1080kg 燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N 的推力，火箭引擎在燃料用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，与火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
来加强对用数值解法解常微分方程实际问题的能力。
关键词：常微分方程数值解 MATLAB目录摘要i绪论1问题解答 (2)总结 (17)参考文献资料 (17)绪论很多科学技术和工程问题常用微分方程的形式建立数学模型，因此微分方程的求解是很有意义的。
建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。
如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些典型的方程，而对于绝大多数的微分方程问题，很难或者根本不可能得到它的解析解，实际问题终归结出来的微分方程主要靠数值解法。
因此，研究微分方程求解的数值方法是非常有意义的。
问题解答问题一：比拟Adams 四阶PECE 模式和PMECME 模式。
数学实验课程设计常微分方程数值解
数学实验报告1.题目：某容器盛满水后，底端直径为d0的小孔开启（如图1），根据水力学知识，当水面高度为h时，谁从小孔中流出的速度为v=0.6*（g*h）^0.5（其中g为重力加速度，0.6问哦小孔收缩系数）1）若容器为倒圆锥形（如图1），现测得容器高和上底直径都为1.2m，小孔直径d为3cm，为水从小孔中流完需要多少时间；2min时水面高度是多少。
2）若容器为倒葫芦形（如图2），现测得容器高1.2m，小孔直径d为3cm，由底端（记x=0）向上每隔0.1m测出容器的直径D（m）如表1所示，问水从小孔中流完需要多少时间；2min时水面高度是多少。
X/m00.10.20.30.40.50.60.70.80.9 1.0 1.1 1.2D/m0.030.050.080.140.190.330.450.680.981.11.21.131.02.分析：由题知，水从小孔中流出，不仅与容器有关，还与水流速度v=0.6*（2*g*h）^0.5有关。
第一小题容器是圆锥形，比较规则，但是由于水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，所以要用到微积分解决。
由（1）知，水面直径等于水深。
水深为h时，流量为0.6(π/4)d^2*(gh)^0.5，0.6*(g*h)^(0.5)*π*(d0/2)^2*dt=π/4*h^2*dh则水深下降dh 所需时间 ：dt=-[(π/4)h^2*dh]/[0.6(π/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]水深由1.2m 至0定积分得水从小孔流完的时间：T （其中已知d=0.03m ，g=9.8m*s(-2）对于第二问：设两分钟（120S ）后水深为X m ，由dt=-[(π/4)h^2*dh]/[0.6*(π/4)*d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]则263.93-120 =X^2.5/[1.5*d^2*(g)^0.5]以d=0.03m ，g=9.8m*s(-2代入上式得 水深：X第二小题容器为倒葫芦形，比较不规则，比较复杂，不仅要考虑水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，所以要用到微积分解决，还要注意表1的倒葫芦形的不断变化，水深的高度变化是不规则的但仍可以用微积分。
实验5常微分方程的数值解
实验5 常微分方程的数值解概要：将装满放射性废物的圆桶扔到水深300ft 的海底，圆桶体积55gal ，装满废料的桶重为527.436lbf ，在海中浮力为470.327lbf 。
此外，下沉时受到的阻力与速度成正比，比例系数为0.08lbf/s 。
实验发现当圆桶速度超过40ft/s 时，就会因与海底冲撞而破裂。
要求：（1）建立解决上述问题的微分方程模型（2）用数值和解析两种方法求解微分方程，并回答谁赢得了官司。
模型建立由牛顿第二定律可列出圆桶下沉速度的微分方程：dv G F f G F bv dt m m ----==其中G 为圆桶重量，F 为浮力，b 为下沉阻力与速度关系的比例系数。
换算到国际单位制，dept=300*0.3048=91.4400 海深（m ）ve=40*0.3048=12.1920 速度极限（超过就会使圆筒碰撞破裂）(m/s) G=527.436*0.4536*9.8＝2344.6 圆筒重量（N ） F=470.327*0.4536*9.8＝2090.7 浮力(N)m=527.436*0.4536＝239.24 圆筒质量（kg ） b=0.08*0.4536*9.8/0.3048＝1.1667 比例系数(Ns/m) 模型求解 一．求数值解Matlab 程序如下： sd.m:function dx=sd(t,x,G,F,m,b) dx=[(G-F-b*x)/m];%微分方程sddraw.m: clear;G=527.436*0.4536*9.8;%圆筒重量（N ） F=470.327*0.4536*9.8;%浮力(N)m=527.436*0.4536;%圆筒质量（kg ）b=0.08*0.4536*9.8/0.3048%比例系数(Ns/m) h=0.1;%所取时间点间隔ts=[0:h:2000];%粗略估计到时间2000 x0=0;%初始条件opt=odeset('reltol',1e-3,'abstol',1e-6);%相对误差1e-6，绝对误差1e-9 [t,x]=ode45(@sd,ts,x0,opt,G ,F,m,b);%使用5级4阶龙格—库塔公式计算 %[t,x]%输出t,x(t),y(t)plot(t,x,'-'),grid%输出v(t)的图形 xlabel('t'); ylabel('v(t)');%用辛普森公式对速度积分求出下沉深度 T=20;%估计20s 以内降到海底for i=0:2:10*T%作图时间间隔为0.2 y=x(1:(i+1)); k=length(y);a1=[y(2:2:k-1)];s1=sum(a1); a2=[y(3:2:k-1)];s2=sum(a2);z4((i+2)/2)=(y(1)+y(k)+4*s1+2*s2)*h/3;%辛普森公式求深度 endi=[0:2:10*T]; figure;de=300.*0.3048.*ones(5*T+1,1);%海深ve=40.*0.3048*[1 1];%速度极限值（超过就会使圆筒碰撞破裂）plot(x(i+1),z4',x(i+1),de,ve,[0 z4(5*T+1)]);%作出速度－深度图线，同时画出海深和速度要求grid;gtext('dept'),gtext('Vmax');xlabel('v');ylabel('dept(v)');figure;plot(i/10,z4');%作出时间－下降深度曲线grid;xlabel('t');ylabel('dept(t)');求解结果如下图：速度—时间曲线：可以看到经过足够长的时间后，如若桶没有落到海底，它的速度会趋于常值。
微分方程数值解法课程试验题目
计算实验课微分方程数值解法数值计算实验题目一、常微分方程部分：1．使用四阶Runge-Kutta 方法求解如下初值问题的近似解，并将结果与实际值进行比较。
2．使用四阶Adams 预估校正算法（PECP 和PMECME 方案），初始值用四阶Runge-Kutta 方法提供，并将结果与实际值进行比较。
()21u t u -+='，32≤≤t ，()12=u ；精度510-=ε，5.0=h 。
实际解11u t t=+-。
tuu +='1，21≤≤t ，()21=u ；精度510-=ε，2.0=h 。
实际解2ln +=t t u 。
二、偏微分方程部分：1．用有限差分法求解如下Poisson 方程(),cos3sin ,u x y x y π-∆=，0x π<<，10<<y边界条件为： ()(),0,10,u x u x ==01x ≤≤； ()()0,,0,x x u y u y π==10≤≤y 取1,h Nπ=和21,h N=作矩形剖分，网格节点为1i x ih =，2j y jh =，i ,j =0,1,…,N 。
差分格式为1,,1,,1,,11222cos3sin i j i j i j i j i j i j i j u u u u u u x y h h π+-+--+-+⎛⎫-+= ⎪⎝⎭,1,2,,1i j N =-L 边界条件为： 00,0,,,i iN u u i N ===L01,1,,1,j j u u j N ==-L 1,1,,1,N j N j u u j N -==-L 结果与精确解()()12,9cos3sin u x y x y ππ-=+进行比较。
求解方案：依次令4,8,16,32N =，取6位小数计算。
用消元法求解，并就(),,44j i j x y π⎛⎫= ⎪⎝⎭，,1,2,3i j =处列出差分解与精确解。
其次，就N =32，0.25，0.5，0.75及i =0，2，4，…，30，32画出差分解曲线。
(整理)微分方程数值解法课程设计报告
《微分方程数值解》课程设计课设题目：团队成员：081110209 丘凯倩081110202 江雨芮 081110232 黄东方 081110310 曲 健 081110311 代永轩指导教师：王春武二〇一四年六月十六日T5-分别利用欧拉公式、改进的欧拉公式和经典的四级四阶龙格-库塔公式求解常微分方程组的初值问题目录一、第十组团队成员及分工 (3)二、研究问题 (4)三、理论分析 (4)四、数值方法 (6)五、计算结果 (9)六、总结及体会 (11)一、第十组团队成员及分工【丘凯倩】081110209 组长统一规划团队课程设计分工，并及时分配任务。
同时，负责给出格式、求解格式的截断误差，对小组成员的成果进行最终汇总。
【江雨芮】081110202 组员分析课设题目，运用MATLAB进行总体编程，负责代码改进与调试工作，并给出各个格式的程序流程图。
【黄东方】081110232 组员分析课设题目，给出“欧拉公式”的求解方法及编程思路，代表小组进行汇报展示。
【曲健】081110310 组员按组长分配的任务，给出“改进欧拉公式”的求解方法及编程思路，并负责PPT及课程设计报告的撰写工作。
【代永轩】081110311 组员配合团队其他成员，给出“四级四阶R-K法”的求解方法及编程思路，并负责PPT及课程设计报告的撰写工作。
二、研究问题分别利用欧拉公式、改进的欧拉公式和经典的四级四阶龙格-库塔公式求解常微分方程组的初值问题：21141(0),2dy y dx y ⎧⎡⎤=⎪⎢⎥-⎪⎣⎦⎨⎛⎫⎪= ⎪⎪⎝⎭⎩01001,.x h <≤=。
三、理论分析1.欧拉公式),(y 1n n n n y x hf y +=+为简化分析，人们常在)(y n n x y =的前提下估计误差*11)(++-n n y x y 。
这种误差称为局部截断误差。
对于欧拉格式，))(,()(*1n n n n x y x hf x y y +=+而按泰勒公式展开有)(2)()()(21ξy hx y h x y x y n n n ''+'+=+因此有)(2)(2*11ξy h yx y n n ''=-++ 所以欧拉公式的截断误差为)(2h O 。
微分方程数值解法数值计算实验题目
y
x
数 值 解 和 精 确 解 的 误 差 曲 面 (N=8)
0.05
u
0
-0.05 1 4 0.5 1 0 3 2 0
y
xቤተ መጻሕፍቲ ባይዱ
数 值 解 和 精 确 解 的 误 差 曲 面 (N=16)
0.02
0.01
u
0
-0.01
-0.02 1 4 0.5 1 0 3 2 0
y
x
数 值 解 和 精 确 解 的 误 差 曲 面 (N=32)
y
x
Δu -0.0006 -0.0000 0.0006 -0.0009 -0.0000 0.0009 -0.0006 -0.0000 0.0006
程序运行中,请稍等…… =================================================================== i 1.0000 1.0000 1.0000 2.0000 2.0000 2.0000 3.0000 3.0000 3.0000 j xi yj u(精确) u1(数值) ------------------------------------------------------------------1.0000 0.7854 0.2500 -0.0265 -0.0271 2.0000 3.0000 1.0000 2.0000 3.0000 1.0000 2.0000 3.0000 0.7854 0.7854 1.5708 1.5708 1.5708 2.3562 2.3562 2.3562 0.5000 0.7500 0.2500 0.5000 0.7500 0.2500 0.5000 0.7500 -0.0000 0.0265 -0.0375 -0.0000 0.0375 -0.0265 -0.0000 0.0265 -0.0000 0.0271 -0.0384 -0.0000 0.0384 -0.0271 -0.0000 0.0271
微分方程数值解法课程设计模板
目录一、问题的描述及算法设计 (1)1.1问题的描述 (1)1.2算法设计 (1)二、算法的理论依据及其推导 (4)2.1二维问题的五点差分格式的理论推导 (4)2.2高斯－赛得尔（Gauss-Seidel）迭代法的理论推导 (5)三、算法的流程图 (6)3.1矩形网格的五点差分流程图 (6)3.2高斯－赛得尔（Gauss-Seidel）迭代法流程图 (7)四、相关的数值结果 (8)五、总结 (9)六、附件 (10)二维椭圆问题的差分格式一、问题的描述和算法设计1.1 问题的描述有限差分方法就是一种数值解法，它的基本思想是先把问题的定义域进行网格剖分，然后在网格点上，按适当的数值微分公式把定解问题中的微商换成差商，从而把原问题离散化为差分格式，进而求出数值解。
此外，还要研究差分格式的解的存在性和唯一性、解的求法、解法的数值稳定性、差分格式的解与原定解问题的真解的误差估计、差分格式的解当网格大小趋于零时是否趋于真解(即收敛性)，等等。
而G －S 迭代法是用逐次逼近的方式得到差分方程组的解，它的存储量小，程序简单，计算量小，因此常用于差分方程组的求解。
偏微分方程边值问题的差分法是物理上的定常问题，其定解问题为各种边值问题， 即要求解在某个区域内满足微分方程，在边界上满足给定的边界条件。
椭圆型方程的差分解法可归结为选取合理的差分网格，建立差分格式，利用G －S 迭代法求解代数方程组。
1.2 算法设计算法Poisson 方程的有限差分方法 求Poisson 方程()()()y x f y x yu y x xu ,,,2222=∂∂+∂∂， dy c b x a ≤≤≤≤,的近似解，其边界条件为()(),,,y x g y x u = 若x ＝a 或x ＝b 且d y c ≤≤ 和()(),,,y x g y x u = 若y ＝c 或y ＝d 且b x a ≤≤输入 端点a 、b 、c 、d ；整数m ≥3，n ≥3；误差要求TOL ；最大迭代次数N 。
实验4_常微分方程数值解
实验4 常微分方程数值解【实验目的】1. 练习数值微分的计算；2. 掌握用MATLAB 软件求微分方程初值问题数值解的方法；3. 通过实例学习用微分方程模型解决简化的实际问题；4. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
【实验内容】题目3小型火箭初始重量为1400kg，其中包括1080kg 燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N 的推力，火箭引擎在燃料用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
odefun=@(t,v)(-((882*t)/5 - (2*v^2)/5 + 18280)/(18*t - 1400));[t,v]=ode45(odefun,[0:0.1:60],[0]);subplot(2,2,1)plot(t,v);grid ontitle('图1.速度-时间')a=diff(v)/0.1;t2 = [0:0.1:59.9];subplot(2,2,2)plot(t2,a);grid ontitle('图2.加速度-时间')s = cumsum(v).*0.1;subplot(2,2,3)plot(t,s);grid ontitle('图3.高度-时间')3.1 燃料燃烧过程物理模型分析设火箭质量为m，高度为h，速度为v，加速度为a，火箭推力为F，重力加速度为g，阻力为f。
1.由火箭上总共携带燃料1080kg，燃料燃烧率为18kg/s，可知火箭上升时间t=60s时，燃料全部烧尽。
2.由阻力正比于速度的平方，比例系数0.4kg/m，可知阻力表达式为f=0.4v2。
3.由于燃料燃烧，火箭的质量是时间的函数，易知m(t)=m0-18t4.5.根据牛顿第二运动定律，有。
【清华】实验4-常微分方程数值解(011813)
实验4 常微分方程数值解化学工程系分9班焦阳2009011813 【实验目的】1. 掌握用MATLAB软件求微分方程初值问题数值解的方法；2. 通过实例学习用微分方程模型解决简化的实际问题；3. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
【实验内容】1.题目3：小型火箭初始重量为1400kg，其中包括1080kg燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃烧用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
建立模型并进行分析：假设火箭在上升过程中，重力加速度g不随高度而变化，即固定g = 9.8m/s^2。
、(1)从火箭开始上升到引擎关闭：设火箭质量为m，高度为h，速度为v，加速度为a，阻力为f：，，ﭸ由牛顿第二定律可得：总ﭸ综上可得：；ﭸ；初值条件为：，；定义域为：。
根据常微分方程组的初值问题，在MATLAB中计算数值解，记，，， 。
通过解出微分方程的数值解，并进行绘图得到高度-时间曲线，速度-时间曲线，加速度-时间曲线如下：由MA度、速度t(s MATLAB 计算速度、加速度如下t(s) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 计算得到的火箭度如下表：h(m)06.57326.4459.76106.5166.7240.2326.7425.7536.9659.8793.6937.81091.1254.1425.1604.1790.1983.2181.2384.的火箭从开始上h(m) 0 6.5737 26.444 59.762 106.57 166.79 240.27 326.72 425.79 536.99 659.8 793.63 937.85 1091.8 1254.7 1425.9 1604.8 1790.8 1983.1 2181.2 2384.5 开始上升到关闭引v(m/s)0 13.18926.57740.06253.53566.89 80.02192.829105.22117.11128.43139.14149.18158.55167.23175.22182.55189.22195.27200.75205.7 关闭引擎这段时间/s) 189 577 062 535 021 829 .22 .11 .43 .14 .18 .55 .23 .22 .55 .22 .27 .75 段时间内各时刻a(m/s^2)13.0571413.304513.4532813.4971913.4331313.2613112.9853412.6121912.1519511.6169311.0212710.38 9.7083329.0209048.3309057.6502496.9900526.3593395.7646125.2094884.694626各时刻的高^2)71404532871931313153421919569312733290490524905233961248862621 2592.4 210.18 4.22201222 2804.5 214.19 3.79432623 3020.6 217.79 3.41201724 3240.1 221.01 3.07303925 3462.7 223.92 2.7726326 3687.9 226.56 2.50440927 3915.6 228.97 2.26774828 4145.6 231.14 2.06332529 4377.8 233.11 1.88975930 4611.9 234.91 1.74334931 4847.7 236.57 1.6177932 5085 238.14 1.50617933 5323.8 239.61 1.40954434 5564.1 240.99 1.32933135 5805.8 242.28 1.2650336 6048.7 243.5 1.21393737 6292.9 244.68 1.17083938 6538.1 245.83 1.1302639 6784.5 246.96 1.09469840 7032 248.05 1.06634841 7280.5 249.1 1.04558942 7530.2 250.12 1.03075943 7780.9 251.14 1.01781844 8032.5 252.15 1.00237545 8285.1 253.16 0.98756846 8538.8 254.15 0.97632347 8793.4 255.12 0.96963948 9049 256.07 0.96629149 9305.6 257.03 0.96239250 9563.1 257.99 0.95272351 9821.5 258.95 0.9411852 10081 259.9 0.93372753 10341 260.83 0.93277654 10603 261.75 0.93633455 10865 262.67 0.9369556 11128 263.61 0.92584757 11392 264.54 0.91379358 11657 265.46 0.91062859 11923 266.35 0.9161260 12190 267.26 0.917011根据表格可以很容易得到：关闭引擎的瞬间，h=12190m，v=267.26m/s，a=0.917011m/s^2。
微分方程数值解问题复习题
dy = λ y 运用这些格式。作为课程设计问题之一，具体的步 dx
骤已经在上课的时候讲过，请自己写上。例如，对于经典四级四阶 Runge-Kutta 格式，我们如此求其绝对稳定区域。 经典四级四阶 Runge-Kutta 格式为
1 ⎧ ⎪ yn +1 = yn + 6 h( K1 + 2 K 2 + 2 K 3 + K 4 ) ⎪ ⎪ K1 = f ( xn , yn ) ⎪ 1 1 ⎪ ⎨ K 2 = f ( xn + h, yn + hK1 ) 2 2 ⎪ 1 1 ⎪ ⎪ K 3 = f ( xn + 2 h, yn + 2 hK 2 ) ⎪ ⎪ ⎩ K 4 = f ( xn + h, yn + hK 3 )
3
⎧ ⎧ ⎪1 − c1 − c2 = 0 ⎪c1 + c2 = 1 ⎪ ⎪ 1 ⎪1 ⎪ 3 根据 en +1 = O(h ) ，必须 ⎨ − a2 c2 = 0 ，也就是 ⎨a2 c2 = 。 2 ⎪2 ⎪ 1 ⎪1 ⎪ − c2b21 = 0 b21c2 = ⎪ ⎪ ⎩2 2 ⎩
1 令 c1 = c2 = , a2 = b21 = 1 ，就得到了预报-校正格式： 2 1 ⎧ ⎪ yn +1 = yn + 2 h( K1 + K 2 ) ⎪ ⎨ K1 = f ( xn , yn ) ⎪ K = f ( x + h, y + hK ) n n 1 ⎪ 2 ⎩ 6.求二级二阶，三级三阶，四级四阶 Runge-Kutta 格式的绝对稳定区域。(分别选
⎛t ⎞ ⎛ t ⎞ t (t − 1) ⋅⋅⋅ (t − j + 1) ⎛t ⎞ ，特别地， ⎜ ⎟ = 1 ， ⎜ ⎟ = t 。 ⎜ ⎟= j! ⎝0⎠ ⎝ 1⎠ ⎝ j⎠ a j = (−1) j ∫ −t (−t − 1) ⋅⋅⋅ (−t − j + 1) 1 1 dt = ∫ t (t + 1) ⋅⋅⋅ (t + j − 1)dt 0 j! j! 0
第六章_常微分方程初值问题的数值解法_习题课
h2 h3 y ( x n ) y ( x n ) O(h 4 ) 2 6 而且 y ( x n ) f ( x n , y ( x n )) ， y ( x n 1 ) f ( x n 1 , y ( x n 1 )) ，对 y ( x n 1 ) 也在 x n 处作 Talor 展开， y ( x n 1 ) y ( x n ) hy ( x n )
湖北民族学院理学院《数值计算方法》教学辅导材料
陈以平编写
h2 h3 y ( x n ) y ( x n ) O(h 4 ) 2 6 h h h2 h3 y ( x n ) y ( x n ) y ( x n ) y ( x n ) y ( x n ) O(h 4 ) 2 2 2 12 h3 y ( x n ) O(h 4 ) O(h 3 ) 12 h3 所以，梯形公式是 2 阶方法，其截断误差的主项是 y ( x n ) 。 12 y ( x n ) hy ( x n )
y k (0.9 0.1y k sin x k ) 0.1( y k 1 y k 1 sin x k 1 )
2
当 k=0，x0=1, y0=1 时，x1=1.2，有 y y (. . y sin x ) (. sin ) .
y f ( x, y ) 3.求解初值问题 欧拉法的局部截断误差是( y ( x ) y 改进欧拉法的局部截断误差是( ); 四阶龙格－库塔法的局部截断误差是( ). (A)O(h2) (B)O(h3) (C)O(h4) (D)O(h5)
4. 改进欧拉法的平均形式公式是( ) y p y k hf ( x k , y k ) y p y k hf ( x k , y k ) (B) y c y k hf ( x k , y p ) .(A) y c y k hf ( x k , y p ) y k ( y p y c ) y k ( y p y c ) y p y k hf ( x k , y k ) y p y k hf ( x k , y k ) (C) y c y k hf ( x k , y p ) (D) y c y k hf ( x k , y p ) y k h ( y p y c ) y k ( y p y c ) (D) 答案：
微分方程数值解 课程设计
微分方程数值解课程设计姓名 *****学号 200******专业信息与计算科学课设题目：对初边值问题2222xu t u ∂∂=∂∂ (0<x<1) 0||10====x x u u 0>tx x u πsin 81)0,(=0|0=∂∂=t t u 10≤≤x 利用显式差分格式：)2(211211n m n m n m n m n m n m U U U r U U U -+-++-=+-求近似解。
其中1,025.0==∆=∆T x t ，并画出图形与解析解比较计算结果。
方程的解析解为)]}(sin[)]({sin[161t x t x u ++-=ππ 算法描述：1. 网格剖分 取]1,0[],1,0[∈∈x t)......2,1,0(，)......2,1,0(，n j jt t t m i ih x x j i ======2. 差分格式)2(211211n m n m n m n m n m n m U U U r U U U -+-++-=+- ...... (1)3. 初值处理,010====x x u u ,0),(),0(==⇒j m u j u ...... (2),0|0=∂∂=t tu ,0)0,()1,(=-⇒i u i u ...... ...... (3)x x u πsin 81)0,(=)),*(sin(*8/1)0,(i x i u ∏=⇒ ……（4） 由（1）~（4）通过迭代逐层可以求得数值解.j i u4．实验结果时间、空间均为0~1，且网格为4141⨯的数值与图像结果： >> x=[0:0.025:1];>> y=[0:.025:1];>> mesh(x,y,u)>> mesh(x,y,u1)近似解图像：精确解图像表一部分数值结果部分数值结果比较图像0.10.20.30.40.50.60.70.80.91Matlab程序：u=zeros(41,41) u1=zeros(41,41)h=0.025% ****³初值**************for i=1:41u(i,1)=0u(i,41)=0endfor j=2:40u(41,j)=(1/8)*[sin(pi*j*h)]u(40,j)=(1/8)*[sin(pi*j*h)]end%************近似解**************for i=39:-1:2for j=2:40u(i,j)=u(i+1,j-1)+u(i+1,j+1)-u(i+2,j)endend%**************精确解***************for i=1:40for j=1:40x1=i*h;t1=j*h;u1(i,j)=1/16*sin(pi*(x1-t1))+1/16*sin(pi*(x1+t1)) endend。
微分方程数值解法(李荣华3版)第二章习题答案(大)
第二章习题课（2007.4.28）习题1．求两点边值问题22sin , 0142(0)0, (1)0xLu u u x u u ππ⎧''=-+=<<⎪⎨⎪'==⎩（1.1）的线性有限元解函数（区间等距剖分成2段或3段），要求在计算总刚度矩阵和总荷载向量时，所涉及的定积分用两种方法： 1. 精确求解；2. 用中矩形公式近似计算。
解：第一步：写出原问题（1.1）的等价变分形式（基于虚功原理）试探函数空间和检验函数空间均为：11(){ |(), ()0 }E H I u u H I u a =∈=.在（1.1）的第一个式子两边同时乘以检验函数空间1()E H I 中的任意元素v ，再在区间(0,1)I =上积分，可得21112sin42xu vdx uvdx vdx ππ''-+=⎰⎰⎰ （1.2）其中111011[(1)(1)(0)(0)]u vdxu v dx vu u v dx v u v u u v dx'''''-=-''''=--''=⎰⎰⎰⎰分部积分（1.3）将（1.3）代入（1.2），可得211()2sin42xu v uv dx vdx ππ''+=⎰⎰记21010(,)()4()2sin 2a u v u v uv dx x f v vdxππ⎧''=+⎪⎪⎨⎪=⎪⎩⎰⎰ 则可以得到原问题（1.1）的等价变分问题：求1()E u H I ∈，使得1(,)(), ()Ea u v f v v H I =∀∈. (1.4)第二步：线性有限元空间的构造1.网格剖分（这里以等距剖分3段为例）2.一次Lagrange 有限元空间的定义1{ ():|(),1,2,3, (0)0 }E i h h h e i h V u C I u P e i u =∈∈==.3. Lagrange 节点基函数的构造113, [0,]312()23, [,]330,x x x x x φ⎧∈⎪⎪⎪=-∈⎨⎪⎪⎪⎩在别处 ； 21231, [,]332()33, [,1]30,x x x x x φ⎧-∈⎪⎪⎪=-∈⎨⎪⎪⎪⎩在别处； 3232, [,1]()30,x x x φ⎧-∈⎪=⎨⎪⎩ 在别处.4.空间E hV 中元素的（整体）表示记 (), 1,2,3i h i u u x i ==，则对E hh u V ∀∈，有31()()h j j j u x u x φ==∑ (1.5)第三步：写出线性有限元方程将原变分问题（1.4）中1()EHI 的试探函数子空间和检验函数子空间均取为E h V ，则可以得到原问题（1.1）的近似变分问题：求 E hhu V ∈，使得 (,)(), E h h h h h a u v f v v V =∀∈. (1.6)利用（1.5）并将 h v 取为(), 1,2,3i x i φ=则上述近似变分问题等价于求123,,u u u R ∈，使得31(,)(), 1,2,3j j i i j a u f i φφφ===∑⇔ 31(,)(), 1,2,3j i j i j a u f i φφφ===∑⇔ 31(,)(), 1,2,3i j j i j a u f i φφφ===∑ 写成矩阵形式AU b =其中111213212223313233(,)(,)(,)(,)(,)(,)(,)(,)(,)a a a A a a a a a a φφφφφφφφφφφφφφφφφφ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦，123u U u u ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦， 123()()()f b f f φφφ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦其中（a ） 精确求解以11(,)a φφ和1()f φ的计算为例：212211110122222223311111031222222233103(,)[()]4[()][()]44[3(3)][(3)(23)]44a dxdx dxx dx x dx πφφφφππφφφφππ'=+''=+++=++-+-=⎰⎰⎰⎰⎰1221(,)(,)a a φφφφ==，1331(,)(,)a a φφφφ==，22(,)a φφ=2332(,)(,)a a φφφφ==，33(,)a φφ=11101233103()2sin2 2sin (3)2sin (23)22xf dxx x x dx x dx πφφππ==+-=⎰⎰⎰（b ）中矩形公式近似求解中矩形公式：()()()2baa bg x dx b a g +≈-⎰.以11(,)a φφ和1()f φ的计算为例：222222112221111(,)[3(3)][(3)(23)]34634211 (9)(9)3163162 (9)316a ππφφπππ≈++-+-=+++=+ 111111162()2sin (3)2sin (23)32632222 sin sin32438f ππφππ≈+-=+习题2．导出下面边值问题1122(), ()(), ()()d du Lu p qu f a x bdx dx u a u a u b u b αβαβ⎧=-+=<<⎪⎨⎪''+=+=⎩ （2.1）的线性有限元方程。
常微分方程数值解法课程设计报告
课程设计报告课程设计题目：常微分方程数值解法学生姓名:专业：班级：指导教师：时间：题目：常微分方程数值解法用欧拉方法、改进欧拉方法、3阶龙格－库塔法以及4阶龙格－库塔法求解常微分方程初值问题：1sin (0)==dy y x y dxe一、摘要在matlab 环境下熟悉的运用计算机编程语言并结合龙格－库塔法的理论基础对常微分方程初值问题进行求解，在运行完程序后以及对运行结果做出各方面的分析和比较。
二、设计目的用熟悉的计算机语言编程上机完成用欧拉方法、改进欧拉方法、3阶龙格－库塔法以及4阶龙格－库塔法求解常微分方程初值问题。
二、理论基础1.欧拉公式【3】在点n x 将1()()n n y x y x h +=+作Taylor 展开，得211()()'()''(),(,)2!n n n n n n n h y x y x hy x y x x ξξ++=++∈，那么当h 充分小时，略去误差项2''()2!n h y ξ，用n y 近似替代()n y x 、1n y +近似替代1()n y x +，并注意到'()(,())n n n y x f x y x =，便得0010(),(,),,0,1,,1,.n n n n n y y x y y hf x y x x nh b a n N h N +⎧⎪=⎪=+=+⎨⎪-⎪=⋅⋅⋅-=⎩上述方法称为Euler 方法。
2.改进Euler 方法【3】在应用梯形方法的迭代公式进行运算时，每迭代一次都要重新计算函数(,)f x y 的值，且还要判断何时可以终止或转下一步计算。
为了控制计算量和简化算法，通常只迭代一次就转入下一步计算。
具体说，我们先用Euler 公式求得一个初步的近似值1n y +，称之为预测值，然后用梯形方法的迭代公式作一次迭代得1n y +，即将1n y +校正一次，这样建立的预测-校正方法称之为改进的Euler 方法： 预测：1(,)n n n n y y hf x y +=+， 校正：111[(,)(,)]2n n n n n n h y y f x y f x y +++=++. 3.三阶龙格-库塔方法【3】类似前面改进的Euler 方法公式的推导方法，将1n y +在(,)n n x y 处作Taylor 展开，然后再将1()n y x +在n x x =处作Taylor 展开，只要将两个展开式前四项相同便有411()()n n y x y o h ++-=。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
课程设计题目
1.计算并画出二级二阶、三级三阶和四级四阶Runge-Kutta 方法的绝对稳定区域。
2.用古典隐式格式并结合追赶法计算抛物型方程
2201,00.01(,0)sin()01
(0,)(1,)00u u x t t x u x x x u t u t t π⎧∂∂= <<<<⎪∂∂⎪⎪= ≤≤⎨⎪== ≥⎪⎪⎩
的近似解。
其中0.1,0.001h k ==。
并与解析解
2sin t u e
x ππ−=
比较。
3.用五点差分格式
,1,1,,1,1,2
1(4)0l m l m l m l m l m l m U U U U U U h +−+−◊=
+++−= 求解Dirichlet 问题 {}2222220(,),(,)0,1ln (1)u u x y x y x y x
y u x y ∂Ω
⎧∂∂+= ∈ΩΩ=<<⎪∂∂⎨⎪⎡⎤=++⎣⎦⎩ 其中，1.0=Δ=Δy x ，分别用Jacobi 迭代和Gauss －Seidel 迭代求解，误差为310−，分析两方法的收敛速度，将计算结果画出图形。
4.方程
2004,012,1(,0)1,1u u x t t x x u x x ∂∂⎧+= <<<<⎪∂∂⎪⎨≤⎧⎪=⎨⎪>⎩⎩
其中0.1h =，0.01k =。
分别用Lax-Friedrichs 格式、C-I-R 格式和Lax-Wendroff 。
格式求解上述微分方程初边值问题。
并画出图形比较计算结果。
5.对初边值问题
22220001
01,011sin ,018|0,010t t x x u u x t t
x u x x u x t
u u t π====⎧∂∂= <<<<⎪∂∂⎪⎪=≤≤⎪⎨⎪∂=≤≤⎪∂⎪⎪== 0≤≤1⎩ 利用显式差分格式：
)
2(211211n
m n m n m n m n m n m U U U r U U U −+−++−=+− 求近似解。
其中0.025h k ==，并画出图形与解析解比较计算结果。
方程的解析解为
[][]{}1
sin ()sin ()16u x t x t ππ=−++。