计算方法实验报告习题3(浙大版)

计算方法实验报告

实验名称： 实验3 梯形法解方程组的初值问题 1 引言

2 实验目的和要求

用欧拉法梯形格式求解方程组初值问题（取h =0.1），所需初值用经典R-K 方法计算。将编程运算的结果与习题七中第5题的计算结果作一比较[1]。

3 算法原理与流程图

（1）原理 对于初值问题

α

=∈=)(],[),,(a y b a x y x f dx

dy

先将[a,b]作n 等分，得到各离散节点),,2,1,0(n i ih a x i =+=。式中，n

a

b h -=

。 设)(x y y =为方程的解，则)(1+i x y 在),(i i y x 处的泰勒展开式忽略高阶无穷小量后得

)1,2,1,0))((,()()(1-=+≈+n i x y x hf x y x y i i i i

)(),(1i i x y x y +分别用近似值i i y y ,1+代入，则得

)1,,2,1,0)(,(1-=+=+n i y x hf y y i i i i

即欧拉格式。用差商近似代替导数也可得到欧拉格式。用向后差商表示两点数值微分公式，得到隐式欧拉格式

)1,,2,1,0)(,(111-=+=+++n i y x hf y y i i i i

两种方法公式相加再取平均值得到梯形格式

[])1,,2,1,0(),(),(2

111-=++=+++n i y x f y x f h

y y i i i i i i

专业：电气工程及其自动化

姓名： 李X

（2）流程图

输入系数矩阵u,边界a,b,步

长h,和初值矩阵v

n←(b-a)/n

y j(1)←v(j)

j=1,2,…,m

A=E m-h/2*u

B=E m+h/2*u

i=1,2,…,n

y i(j+1)=A-1B*y i(j)

j=1,2,…,m

i

输出x,y i

4程序代码及注释

%欧拉法梯形格式

function [x,y]=euler(u,v,a,b,h)

%说明

%对于下列类型的一阶方程组进行求解

%y1'=u11*y1+u12*y2+…+u1m*ym;

%y2'=u21*y1+u22*y2+…+u2m*ym;

%…

%ym'=um1*y1+um2*y2+…+umm*ym;

%可写成矩阵乘积的形式，u=[u11,u12,…,u1m;u21,u22,…,u2m;……];

%初值向量v

%y1(1)=v1,y2(1)=v2,…,ym(1)=vm;

n=(b-a)/h;m=length(u);

x=[a:h:b];

y(:,1)=v;

%定义y(:,1)代表方程组解矩阵的初值

p=eye(m);

A=-h/2*u+p; B=h/2*u+p;

%由梯形法公式，推导出m 元一阶方程组的解法 for j=1:n

y(:,j+1)=A\B*y(:,j); end

%解矩阵y 的第j 列为x=a+(j-1)*h 时的数值解 x=x';

%输出列向量 y=y';

%转置后，第i 行第j 列为x=a+(i-1)*h 时yj 的数值解

注释

对于一般形式的m 元一阶微分方程组

⎥

⎥

⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢

⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡m mm m m m m m m

m m m y y y u u u u u u u u u y y y x f y y y x f y y y x f y y y 2121222211121121212211''2'1),,,,(),,,,(),,,,(，满足初值条件⎥⎥⎥⎥⎦

⎤

⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡m m v v v y y y 2121)1()1()1(，

[]b a x ,∈

可以用梯形法求解

⎥

⎥⎥

⎥⎦

⎤

⎢⎢⎢⎢⎣⎡+++⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎦

⎤

⎢⎢⎢⎢⎢⎢⎢⎣⎡+++=⎥

⎥⎥

⎥⎦⎤

⎢⎢⎢⎢⎣⎡++++

++++

++++

+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++)1()1()1(2)()()(122221222212))1(,),1(),1(,())(,),(),(,())1(,),1(),1(,())(,),(),(,())1(,),1(),1(,())(,),(),(,(2)()()()1()1()1(212

1

22221112112121

22221

112112*********

2112112121i y i y i y u u u u u u u u u h i y i y i y u h u h u h u h u h u h u h u h u h i y i y i y x f i y i y i y x f i y i y i y x f i y i y i y x f i y i y i y x f i y i y i y x f h i y i y i y i y i y i y m mm m m m m m mm m m m m m m m m m m m m m m

整理后得到

⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎦

⎤

⎢⎢⎢⎢⎢⎢⎢⎣⎡+++=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡---)()()(12

2221222212)1()1()1(2122

2212

22

21212

122221

112

112121

222

2111211

i y i y i y u h u h u h

u h u h u h u h

u h u h

i y i y i y u h u h u h u h u h u h u h u h u h

m mm m m m m m mm m m m m

记为