4、问题三的求解过程
4.1 问题三的分析
必须利用matlab内部函数dsolve求解此常微分方程,此处取a=-40。 4.2 问题三的求解程序(见附录) 4.3 问题三的求解结果
6x 102154321000.10.20.30.40.50.60.70.80.91
B题的实验过程
1、实验中问题的重述
?du?dt??2000u?999.75v?1000.25??dv?u?v求解微分方程组?任选一显式方法,编写程序取不同步长求?dt?u(0)?0,v(0)??2??解上述方程。
2、实验分析
分析该题题意可得到,该题是一个已知初值与未知量所满足的两个微分方程所组成的一阶微分方程组。需要我们做的就是去应用或构造一种Euler显式方法对其微分方程组进行求解。
分析该一阶微分方程组的求解方法,引进tn?t0?nh, n=1,2,···,以un和vn表示节点tn上的
近似值,此处应用Euler方法中的梯形格式形式进行求解,即将显式Euler格式与隐式Euler格式的
结合,是其两者的算术平均。这样既在一定程度上提高了计算精度,也缩减了计算过程的复杂性,计算简单可行。
3、实验求解过程
首先,确定初值。由已知已经得到u(0)=0;v(0)=-2。除此之外,对向量t,u,v赋初值。假定所需得到的结果个数为n,则将t,u,v均赋初值为一个含有n+1个元素均为0的向量。
然后,确定步长。假定所求方程组的值的范围是在[a,b]内,所需得到的结果个数为n,则可得到步长h为(b-a)/n。
其次,引用梯形格式。令f(u,v)=-2000u+999.75v+1000.25;g(u,v)=u-v。引进初值问题:
?dy??f(x,y),其中n=1,2,···,以un和vn表示节点tn上的近似值,则其梯形公式具有以下?dx??y(x0)?y0形式:
huk?1?uk?[f(un,vn)?f(un?1,vn?1)];
2hvk?1?vk?[g(un,vn)?g(un?1,vn?1)]
24、实验求解结果
通过对梯形公式的运用,在Matlab的命令窗口输入a,b,n的具体取值,则可运算出u,v的具体数值值,即可得到向量u和v,见下表:
表1 参数值
t 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 u 0 0.05 0.22371875 -11.54774764 1129.99566 -109081.7181 10531931.71 -1016864424 98178882731 -9.47923E+12 9.15226E+14 v -2 -1.8975 -1.788939063 -2.265693554 53.76998672 -5346.504634 516063.3204 -49826364.48 4810765869 -4.64482E+11 4.48461E+13 C题的实验过程
1、实验中问题的重述
?dy22x??y?xe,1?x?22x给定初值问题?dxx其精确解为y?x(e?e)。试用改进的欧拉法,步
??y(0)?1长h?0.05,h?0.1,编写程序求在节点xk?1?0.1k(k?1,2,?10)处的数值解及误差,并用MATLAB中的内部函数dsolve求此常微分方程初值问题的解并与上述结果进行比较。
2、实验分析
该问题给定初值求解常微分问题,并用改进的欧拉法计算,且需利用不同的步长进行处理比较,其主要目的是研究步长对求解常微分问题的影响程度.
3、实验求解过程(代码见附录)
利用MATLAB编程实现求解过程,程序见附录。
4、实验求解结果
步长和数值解的比较如下:
x y 1.5 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.55 0.15 0.34 0.58 0.86 1.20 1.60 2.07 2.61 3.24 1.6 1.65 1.7 1.75 1.8 1.85 1.9 1.95 2 3.96 4.78 5.71 6.76 7.95 9.28 10.77 12.44 14.30 16.36 18.65
x y 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 0.34 0.86 1.59 2.60 3.94 5.68 7.91 10.72 14.24 18.58
步长不同对误差的影响: 单位(E-03)
x y 1.5 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 3.003 2.905 2.808 2.714 1.55 1.6 1.65 1.7 2.622 1.75 2.532 1.8 2.445 1.85 2.360 1.9 2.279 1.95 2.200 2.124 2.050 1.979 1.910 1.844 1.779 1.718 1.658 1.600
x y 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 10.24 9.609 9.000 8.419 7.869 7.351 6.864 6.407 5.979
由以上数据可知步长越小其所对应的误差越小。
D题的实验过程
1、实验中问题的重述
?dy22x??y?xe,1?x?22x给定初值问题?dxx其精确解为y?x(e?e)。试用欧拉法,步长
??y(0)?1h?0.025,h?0.1,编写程序求在节点xk?1?0.1k(k?1,2,?10)处的数值解及误差,并用
MATLAB中的内部函数dsolve求此常微分方程初值问题的解并与上述结果进行比较。 2、 实验分析
这道题主要考察用欧拉法编写程序求此常微分方程初值问题,对于误差分析,主要采用局部截断误差,因为实验中的累积误差较复杂。方程中含有导数项,这是微分方程的本质特征,导数是极限过程的结果,而计算过程则总是有限的,所以数值解法的第一步就是消除导数项dy/dx,即y’,这个过程也叫离散化。
?dy??f(x,y)算法原理:对于初值问题?dx假设第n步在点xn的值计算没有误差,即yn?y(xn),
??y(x0)?y0由单步法计算出yn?1,则 Tn?1?y(xn?1)?yn?1称为点xn?1上的局部截断误差。
欧拉法:若用差商
y(xn?1)?y(xn)替代其中的导数项y’(xn),则有近似关系式:
hy(xn?1)?y(xn)?hf(xn,y(xn))
若用yn替代y(xn),这样设计出的计算公式为:
yn?1=yn+hf(xn,yn)
则该初值问题变为:
yn?1= yn+hf(xn,yn); y0?y(x0)
局部截断误差为:
y(x?h)?y(x)?h?(x,y(x),h)??121hy??(x)?h3y???(x)???O(h2) 26所以欧拉方法是一阶方法,其局部截断误差为O(h2)。 3、实验求解过程
分别令h=0.025和h=0.01,带入欧拉公式,当h=0.025和h=0.01时,求出x=(1:h:2) 对应的yn的值,并画出图像。
打开matlab软件,新建一个空白M文件,输入如下程序并保存(见附录);然后,在command window
中输入以下程序,并画出图像如下:
函数真实值与欧拉步长为h=0.025和h=0.01的误差:
4、实验求解结果
结果分析:由以上结果可知:采取欧拉法,步长取得越小,误差就越小;步长越大,误差也越大。但总的分析它的误差值,欧拉法的精度是比较低的,误差相对也较大,这种方法是不适用的。 若用MATLAB中的内部函数dsolve求此常微分方程初值问题的解:打开matlab软件,输入desolve函数,发现其中运行错误。
附 录
A题的求解程序 %%4阶经典R-K法 function R = Rungkuta4(f,a,b,N,ya) % a,b 是自变量的取值区间[a,b]的端点 % N是区间等分的个数 % ya表示初值y(a) % E = [x',y']是自变量x和解y所组成的矩阵 h = (b-a)/N; x = zeros(1,N+1); y = zeros(1,N+1); x = a:h:b; y(1) = ya; for i = 1:N k1 = feval(f,x(i),y(i)); k2 = feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3 = feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4 = feval(f,x(i)+h,y(i)+h*k3); y(i+1) = y(i)+(h/6)*(k1+2*k2+2*k3+k4); end R = [x',y']; (问题1)求解程序(Matlab10.0软件) a= 0; b= 1; a= 0; b= 1;