数值分析课程课程设计

课程设计

我再也回不到大二了，

大学是那么短暂

设计题目数值分析

学生姓名李飞吾

字亏xxxxxxxx

专业班级信息计xxxxx班

指导教师_________________________

设计目的：

通过不同题目的理解，进行算法分析。 通过MATLAB 软件进行编程对题目进行解决。 个别题目设计验证，加深对数值分析的理解。 函数的图像绘制的运用 设计题目： ① 题1. 1 利用逆向递推的方法求解问题，通过条件终止地推 ② 题1. 2

从某个初始值开始，利用递推公式进行积分估值 @题1. 3 绘制Koch 分形曲线，节点之间的关系与坐标变换

(4题2. 1

用高斯消元法的消元过程作矩阵分解， LU 分解

⑤题2. 2 矩阵分解方法求上题中 A 的逆矩阵，针对不同的 b,而重复利用已知的 LU ?题2. 3 验证希尔伯特矩阵的病态性 ，矩阵基本运算

⑦ 题3. 1 用泰勒级数的有限项逼近正弦函数，由图像观察逼近效果 ⑧ 题3. 2 绘制飞机的降落曲线，线性方程组求解，与绘图 ⑤题4. 1 线性拟合的函数表达式的推导，使用了两种代码方法 面题5. 1 用几种不同的方法求积分，观察数值积分的逼近效果 ?题5. 5

求空间曲线L 弧长。求导后使用符号函数积分计算

O 题6. 1 用欧拉公式和四阶龙格-库塔法分别求解下列初值问题，代码搜索内容。 (0题6. 4

常微分方程的解，dsolve()函数使用 @)题8. 2

差分法解常微分方程边值问题，

ode 函数无能为力，Matlab 中提供bvp

食军算器。 solinit=bvpinit(x,yinit,params)sol= bvpsolver(odefun,bcfun,solinit,options) (g 题8. 3 求解圆的半径，圆心。线性方程组解参数

设计总结：

(1) 算法是题目的解题核心，好的算法可以使计算更加精确 (题5.1)

(2) 图形绘制在今后的课程设计，或者是论文中可以用到。 (3) 无法解决的问题，需要请教室友，或者上网查阅。

(4)

MATLAB 是一个很强大的软件，提供了很多内置的数学函数，直接进行解题。 查阅资料时了解到一些 MATLAB 论坛。通过帖子阅读，了解到了 MATLAB 在科

学计算方面，模拟，图形，视频等起到的作用。增加了对“计算科学“的理解。

建议：从学生的工作态度、工作量、设计(论文的)创造性、学术性、使用性及书面表达能力等方面给出评价。

签名： 20 年 月 日

课 程 设 计 主 要 内 容

指 导 老 师 评 语

数值分析谋程设计

1. 1水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的

一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？ ( 15621)试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题解：算法

分析：解该问题主要使用递推算法，关于椰子数目的变化规律可以设起初的椰子数为p o,第一至五次猴子在夜里藏椰子后，椰子的数目分别为p o,p1,p2,P3,P4再设最后每个人分得x个椰子，由题：

R 1 7( p k 1)(k=0,1,2,3,4 ) x -(p5 1)

5 5

所以P5 5x 1, p< p< 1 1利用逆向递推方法求解

R 5妇1, (k=0,1,2,3,4 )

4

MATLAB码：n=input('n='); n= 15621 Array

for x=1:n p=5*x+1; for k=1:5 p=5*p/4+1; end

if p==fix(p), break end end

disp([x,p])

d n

1 x

1. 2 设，I n----------------- dx

05 x

(1) 从10尽可能精确的近似值出发，利用递推公式：

1 I n 5I n 1 —(n 1,2,L 20) n

计算机从I I到I 20的近似值；

(2) 从I30较粗糙的估计值出发，用递推公式：

I n i -I n — (n 30,29,L ,3,2)

5 5n

计算从I l到I20的近似值；

解：首先第一步，估计I o和I30的值：

syms x n;

int (x A0/(5+x),0,1)

ans=log(2)+log(3)-log(5)

eval(ans)

ans=

0.1823

则取I。为0.18

syms x n;

int(xA30/(5+x),0,1) ans =

931322574615478515625*log(2)+931322574615478515625*log(3)-931 322574615478515625*log(5)-79095966183067699902965545527073/46 5817912560

eval(ans)

ans =

MATLAB中小数点后保留四位，由上面计算知道积分的值不为了零。

所以I30的取值为0.00001-0.0001

MATLAB码：

i=input('i=');

i=0.18;

>> if i>=0.1&&i<=0.2

for n=1:1:20

i=-5*i+1/n

end

elseif i>0&&i<=0.0001

for n=30:-1:2

i=(-1/5)*i+1/(5*n)

end

end

i =

1.1336e+005

i =

-5.6679e+005

i =

2.8339e+006

i =

-1.4170e+007

i =

7.0848e+007

i =

-3.5424e+008

i =

1.7712e+009

i =

-8.8560e+009

i =

4.4280e+010

i =

-2.2140e+011

同理输入积分初始值i=0时可以得i=0.0884

结果分析：第二种方法所得的结果相对来说比较精确一些，也比较可靠因为第一种方法每一迭代都将最初的误差放大了五倍，使得最终的误差

越来越大；而第一种方法经过每一次迭代都将误差缩小为初始误差的五分之一，使得最终的误差越来越小，因此相对来说比较可靠，性能较好。

1. 3绘制Koch分形曲线

问题描述：从一条直线段开始，将线段中间的三分之一部分用一个等边

三角形的另两条边代替，形成具有5个结点C0SJ的国写（图1-4）;在新的图形中，又将图中每一直线段中间的询之丁部分都用一个等边三角形的另两条边代替，再次形成新的图形（图si t5）,cos时，图形中共有17

个结点。这种迭代继续进行下去可以形成Koch分形曲线。

问题分析：考虑由直线段（2个点）产生第一个图形（5个点）的过程，设p1和R分别为原始直线段的两个端点。现在需要在直线段的中间依次插入三个点

P2，R，P4产生第一次迭代的图形（图1-4）。显然，只位于P点右端直线段的三分之一处，已点绕p2旋转60度（逆时针方向）而得到的，故可以处理为向量巳已经正交变换而得到向量P2P3,形成算法如下：

(1) P2 P (任P)/3.

(2) R P 2(P5 R)/3.

(3) P3 % (E E) A T.

在算法的第三步中，A为正交矩阵。

这一算法将根据初始数据（P和R点的坐标），产生图1-4中5个结点的

坐标。这5个结点的坐标数组，组成一个 5X 2矩阵。这一矩阵的第一 行为为

R 的坐标，第二行为R 的坐标，第二行为P2的坐标……第五行为P

5 的坐标。矩阵的第一列元素分别为 5个结点的x 坐标，第二列元素分 别为5个结点的y 坐标。

问题思考与实验：

(1) 考虑在Koch 分形曲线的形成过程中结点数目的变化规律。设第 k 次迭代

产生结点数为nk ,第k 1迭代产生结点数为nk 1,试写出nk 和nk 1之 间的递推关系式；

(2) 参考问题分析中的算法，考虑图1-4到图1-5的过程，即由第一次 迭代的5个结点的结点坐标数组，产生第二次迭代的 17个结点的结点 坐标数组的算

法；

(3) 考虑由第k 次迭代的n k 个结点的结点坐标数组，产生第k 1次迭代 的n k 1个结点的结点坐标数组的算法；

(4) 设计算法用计算机绘制出如下的

解：(1) n k1 4n k 3

(2)(3辩法及(4)代码分析：

p=[0 0;10 0]; %P 为初始两个点的坐标，第一列为x 坐标，第二列为y 坐标 n=2; %n 为结点数

A=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)]; % 旋转矩阵

计算相邻两个点的坐标之差，得到相邻两点确定的向量

则d 就计算出每个向量长度的三分之一 ，与题中将线段三等分对应

迭代公式

以原点为起点，前n-1个点的坐标为终点形成向量 迭代后处于4k+1位置上的点的坐标为迭代前的相应坐标 用向量方法计算迭代后处于 用向量方法计算迭代后处于 用向量方法计算迭代后处于 迭代后新的结点数目

2. 1用高斯消元法的消元过程作矩阵分解。设

20 2 3 A 1 8 1

2

3 15

消元过程可将矩阵 A 化为上三角矩阵U,试求出消元过程所用的乘数 m21

、m 31、伉1并以如下格式构造下三角矩阵L 和上三角矩阵U

1 20

2 3

L

m 21 1

,U

a 22) a 23)

4

(2)

Koch 分形曲线(图1-6)

for k=1:4

d=diff(p)/3; %diff % m=4*n-3; % q=p(1:n-1,:); % p(5:4:m,:)=p(2:n,:); % p(2:4:m,:)=q+d; % p(3:4:m,:)=q+d+d*A'; % p(4:4:m,:)=q+2*d; % n=m; % end

plot(p(:,1),p(:,2)) % axis([0 10 0 3.5])

4k+2位置上的点的坐标 4k+3位置上的点的坐标 4k 位置上的点的坐标

m31 m32 1 a33

验证：矩阵A可以分解为L和U的乘积，即A=LU。

矩阵LU分解MATLAB代码：

function hl=zhjLU(A)

[n,n]=size(A);RA=rank(A);

if RA~=n

disp('因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的

秩RA如下：’)；

RA,hl=det(A);

return

end

if RA==n

for p=1：n

h(p)=det(A(1：p,1：p));

end

hl=h(1：n);

for i=1：n

if h(1,i)==0

disp(' 因为A的各阶主子式等不等于零，所以A能进彳f LU

分解.A的秩RA和各阶顺序主子式值hl庆次如下：’)；

RA,hl

return

end

end

if h(1,i)~=0

disp(' 因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA 和各阶顺序主子式值hl如下：’)；

for j=1：n

U(1,j)=A(1,j);

end

for k=2:n

for i=2:n

for j=2:n

L(1,1)=1;L(i,i)=1; if i>j

L(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);

L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k))/U(k,k);

else

U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);

end

end

end

end

RA,hl,U,L

end

end

仅上上代码保存为M文件，并在命令窗口输入

A=[20 2 3;1 8 1; 2 -3 15];

b=[0 0 0]';

h=zhjLU(A)

程序运行结果:

L = U=

1.0000 0 0 20.0000

2.0000

3.0000

0.0500 1.0000 0 0 7.9000 0.8500

0.1000 -0.4051 1.0000 0 0 15.0443

实验验证：可以直接使用MATLAB内置LU分解

A=[20 2 3;1 8 1; 2 -3 15];

[L U]=lu(A) 输出结果与上程序输出结果一致。

2. 2用矩阵分解方法求上题中A的逆矩阵。记

1 0 0

b 0 ,b2 1 ,b3 0

0 0 1分别求解方程组AX b

1，AX b2，AX

由于三个方程组系数矩阵相同，可以将分解后的矩阵重复使用。对第一

个方程组，由于A=LU ，所以先求解下三角方程组LY

知再求解上三角 方程组UX Y,则可得逆矩阵的第一列列向量；类似可解第二、第三方 程组，得逆矩阵的第

二列列向量的第三列列向量。由三个列向量拼装可 得逆矩阵A 1

。 解:MATLAB 代码如下:

b1=[1;0;0]; b2=[0;1;0]; b3=[0;1;1]; A=[20,2,3;1,8,1;2,-3,15];

L=[1,0,0;0.05,1,0;0.1,-0.4051,1]; U=[20 2 3;0 7.9 0.85;0 0 15.0443]; Y1=L\b1 X1=U\Y1

Y2=L\b2 X2=U\Y2 Y3=L\b3 X3=U\Y3 Y1 = 1.0000 -0.0500 -0.1203 X1 = 0.0517 -0.0055 -0.0080 Y2 =

0 1.0000 0.4051 X2 =

-0.0164

0.1237

0.0269

Y3 = 0 1.0000 1.4051

X3 = -0.0257 0.1165 0.0934

实验验证： [X1 X2 X3] ans = 0.0517 -0.0164 -0.0257 -0.0055 0.1237 0.1165 -0.0080 0.0269 0.0934 而：inv (A) =[X1 X2 X3]

得证

2. 3验证希尔伯特矩阵的病态性：对于三阶矩阵

1 1/

2 1/

3 H

1/2 1/3 1/4 1/3

1/4 1/5

取右端向量b 【11/6 13/12 47/6°]T ,验证：

(1) 向量

X [X 1 X 2 X 3

]T [1 1 1]T

是方程组HX b 的准确解；

(2) 取右端向量b 的三位有效数字得 b [1.83

L08 0.783]T

,求方程组 的准确解X ,

并与X 的数据I

1,1,1

】'作比较。说明矩阵的病态性。

解：⑴

H=[1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5]; X=[1；1；1]； b=H*X b =

1.8333 1.0833 0.7833

与题中相同

(2)先求出解X,与数据[〔，"亍作比较 H=[1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5]; b=[1.83;1.08;0.783]; X=H\b X =

1.0800 0.5400 1.4400

与[1,1,1]相差较大，矩阵为病态矩阵

3. 1用泰勒级数的有限项逼近正弦函数

sin x,x [0,] x,x [0, /2] 3

—

—

x x /6,x [0, /2] x x 3/6 x 5/120, x

y °(x ) y 〔(x ) y 2(x ) [0,

用计算机绘出上面四个函数的图形。

解:

MATLAB代码如下(1) syms x; taylor(sin(x)) x=0:0.01*pi:pi plot(x,sin(x))

⑵

syms x;

taylor(x)

x=0:0.01*pi:pi/2

plot(x)

⑶

syms x;

taylor(x-x A3/6)

fplot('x-x A3/6',[0 pi/2]) (4)

syms x;

taylor(x-xA3/6+xA5/120)

fplot('x-xA3/6+xA5/120',[0 pi/2])

结果图形右：

x=0:0.1:pi;

y=sin(x);

plot(x,y,'-sk');

hold on

x=0:0.1:pi/2;

y=x;

plot(x,y,'-b*')

hold on

fplot('x-x.A3/6',[0,pi/2,0,2],2e-3,'-gx')

hold on fplot('x-x.A3/6+x.A5/120',[0,pi/2,0,2],2e-3,'-r.') hold off

legend('sin(x)','x','x-xA3/6','x-A3/6+xA5/120',2)

xlabel('x')

ylabel('y')

title('Taylor approximation')

结果分析：图中红色点线为正弦曲线，蓝色的星线为一阶泰勒逼近，

绿色叉线为二阶

豪勒逼近，黑色正方形线为三阶泰勒逼近。可见三阶泰勒逼近效果最好，泰勒级改越高，逼近效果尬婷。

3. 2绘制飞机的降落曲线

一架飞机飞临北京国际机场上空时，其水平速度为540km/h,飞行高度为1 000m。飞机从距机场指挥塔的横向距离12 000m处开始降落。

根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。建立直角坐标系，设飞机着陆点为原点O，降落的飞机为动点P(x,y)，则x表示飞机距指挥塔的距离，y表示飞机的飞行高度，降落曲线为

2 3

y(x) a0 a1x a2x a3x

该函数满足条件：

y(0) 0, y(12 000) 1000

y'(0) 0, y'(12000) 0

(1) 试利用y(x)满足的条件确定三次多项式中的四个系数；

(2) 用所求出的三次多项式函数绘制出飞机降落曲线。

function S=f(x1,y1,x2,y2,x3,y3,x4,y4) format long

a1=[1 ,x1,x1A2,x1A3];

a2=[1 ,x2,x2A2,x2A3];

a3=[0,1,2*x3,3*x3A2];

a4=[0,1,2*x4,3*x4A2];

a=[a1; a2;a3;a4];

b=[y1; y2;y3;y4];

s=a\b;

x=-12000:250:0;

y=s(3)*x.A2-s(4)*x.A3

plot(x,y,'-d')

xlabel('x')

ylabel('y')

以上为M文件内容，在命令窗口键入

f(0,0,12000,1000,0,0,12000,0)

运行结桌：ans为多项土索数

ans =

1.0e-004 *

0.20833333333333 -0.00001157407407

4. 1曾任英特尔公司董事长的摩尔先生早在 1965年时，就观察到 一件很

有趣的现象：集成电路上可容纳的零件数量，每隔一年半左右 就会增长一

倍，性能也提升一倍。因而发表论文，提出了大大有名的

MATLAB ：码：

x=[1959,1962,1963,1964,1965]; y=[1,3,4,5,6]; p1=polyfit(x,y,1) y1=polyval(p1,x) plot(x,y1,'-',x,y,'r*') xlabel('x'),ylabel('y'); 运行结果： p1 =

1.0e+003 *

0.0008 -1.6255 y1 =

0.8113 3.3019 4.1321 4.9623 5.7925 线性拟合的函数表达式： Y=0.8302x-1.6255e+003 解法二：

x=[1959 1962 1963 1964 1965]; y=[1;3;4;5;6]; for i=1:length(x)

for j=1:2

A(i,j)=x(i)”1); end end

[L,U]=lu(A'*A); xs=U\(L\(A'*y))

从而年代y 与增加倍数x 之间的关系为： y=-1625.528301891238+0.830188679248x

摩尔定律（Moore' s LaW,并预测未来这种增长仍会延续下去。下面 数据中，第二行数据为晶片上晶体数目在不同年代与 1959年时数目

比较的倍数。这些数据是推出摩尔定律的依据： 试从表中数据出发，推导线性拟合的函数表达式。 解：解法一

运行结果： xs =

1.0e+003 *

-1.625528301891238 0.000830188679248

(1)牛顿-莱布尼茨公式；(2)梯形公式；(3)辛卜生公式；(4)复 合梯形公式。

解： syms x

i1=int(4/(1+x A 2),x,0,1) a=0;b=1; h=b-a;

i2=(4/(1+aA2)+4/(1+bA2))/2

i3=h/6*(4/(1+aA2)+4*4/(1+(a+b/2)A2)+4/(1+bA2)) M=100; h=(b-a)/M; i4=0;

for k=1:(M-1)

x=a+h*k;

i4=i4+4/(1+x"2); end

i4=h*(4/(1+aA2)+4/(1+bA2))/2+h*i4

牛顿莱布尼兹公式得到精确结果a 采用挣形公式得到的结果比采用 Simpson 公式的精蒲I 度鱼很，米用复祝榜形公式在步长取得越束裁小 的状态下可以提高精度。

x cost

5. 1用几种不同的方法求积分

―火的值。

运行结果：

11 =

pi

12 =

3 13 = 3.1333 1

4 = 3.1416

5. 5求空间曲线L:

y sint

z 2 cost sint

弧长公式为

L : .x'2

(t)

y'2(t) z'2

(t)dt

解：运行程序为：

syms t;

x=diff(cos(t)); y=diff(sin(t));

z=diff(2-cos(t)-sin(t));

y=int((xA2+yA2+zA2)A0.5,t,0,2*pi); % digits(14); i=vpa(y)

运行结果：i =

8.7377525709894

符号函数积分

6. 1用欧拉公式和四阶龙格-库塔法分别求解下列初值问题;

0.9y

(i)y' , y(0) i;x [0,1]

1 2x

⑵y' 闩顶。)

1 x

解：(1) <1>欧拉公式：

function [t,x]=Euler(fun,t0,tt,x0,N) h=(tt-t0)/N;

t=t0+[0:N]'*h;

x(1,:)=x0';

for k=1:N

f=feval(fun,t(k),x(k,:));

f=f;

x(k+1,:)=x(k,:)+h*f;

end

以'Euler.m'保存

function f=Euler_fun(t,x) f=0.9*x./(1+2*t)「以'Euler_fun.m'保存function

main_Euler

[t,x]=Euler('Euier_fun',0,1,1,20);

fh=dsolve('Dx=0.9*x/(1+2*t)','x(0)=1'); for k=1:8

ft(k)=t(k);

fx(k)=subs(fh,ft(k));

end

[t,x] , 4

以'main_Euler.m'保存输入：main_Euler <2>四阶龙杯-库塔法

function R=rk4(f,a,b,ya,N)

%y'=f(x,y)

%a,b为左右端点

%Nfe迭代步长

%h为步长

%ya为初值

h=(b-a)/N;

T=zeros(1,N+1);

Y=zeros(1,N+1);

T=a:h:b;

Y(1)=ya;

for j=1:N

k1=h*feval(f,T(j),Y(j));

k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);

k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);

k4=h*feval(f,Y(j)+h,Y(j)+k3);

2；x [0,1]

1.0000

1.0450

1.0878

1.1285

1.1676

1.2051

1.2413

1.2762

1.3100

1.3427

1.3745

1.4055

1.4356

1.4649

1.4936

1.5216

1.5490

1.5758

1.6021

1.6278

1.6531

Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;

end

ans=[T' Y']

以'rk4.m'保存

function z=f(x,y)

z=0.9*y./(1+2*x);

以'f.m'保存

输入：rk4(T,0,1,1,8)

6. 4列出函数f(x) [1 ln (x)] (x)在区间[0, e]上的函数值表并作出

它的图象。其中，(x)是初值问题

d ln[ (x)] 2x dx

(0) 0

的解

解

v=dsolve('Dv*log(x)=2*x','v(0)=0','x');

f=(1-log(v))*v

f =

-2*(1-log(-2*Ei(1,-2*log(x))))*Ei(1,-2*log(x)) ezplot(T)

输出结果：

f =

-2*(1-log(-2*Ei(1,-2*log(x))))*Ei(1,-2*log(x))

7.4 一个10次项式

10 9 8 .

P(x) x a1x a2x L

的系数为

[1 a1 a2 …a9 a10]=[1 -55 1320 T81 50 157 773 电02 055 341 693 0

-840 950 0 127 535 76 -106 286 40 632 880 0]

试用多项式的求根指令roots求出该10次方程的10个根，然后修改9次项的系数-55为-56,得新的10次方程，求解新的方程，观

察根的变化是否很显著。

解：

p=[1 -55 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800];

>> roots(p)

ans =

10.6051 + 1.0127i

10.6051 - 1.0127i

8.5850 + 2.7898i

8.5850 - 2.7898i

5.5000 + 3.5058i

5.5000 - 3.5058i

2.4150 + 2.7898i

2.4150 - 2.7898i

0.3949 + 1.0127i

0.3949 - 1.0127i

p=[1 -56 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800];

>> roots(p)

ans =

21.7335

7.3501 + 7.7973i

7.3501 - 7.7973i

4.3589 + 3.3285i

4.3589 - 3.3285i

5.1831

2.4378 + 2.7974i

2.4378 - 2.7974i

0.3949 + 1.0127i

0.3949 - 1.0127i

结果分析：改变系数之后，根的变化显著

8. 2用差分法解常微分方程边值问题：

y y 1 0,x (0,1) y(0) 0,y(1) 0

取h=0.1,xj=jh(j=0,1,2,』,10)求y1, y2, ? y9并与该问题的准确解

y(x) ―e e x e x 1 1 e 1 e

比较，列出各节点处的近似解、准确解和误差。

解析：ode**函数无能为力Matlab中，提供了bvp解算器。

solinit = bvpinit(x, yinit, params) sol bvpsolver(odefun,bcfun,solinit,options)

解题代码：

clear;close;

sinit=bvpinit(0:1,[0,0]);

odefun=inline('[y(2);y(1)-1]','t','y');

bcfun=inline('[ya(1);yb(1)]','ya','yb sol=bvp4c(odefun,bcfun,sinit)

t=linspace(0,1,101);

y=deval(sol,t);

plot(t,y(1,:),sol.x,sol.y(1,:),'o') legend('解曲线','解点',2) xlabel('x'),ylabel('y') 运行结果：

sol =

x: [0 0.2500 0.5000 0.7500 1] y: [2x5 double] yp: [2x5 double] solver: 'bvp4c'

8. 3通过平面上的三个点P1 (x1,y1)尸2(x2,y2),P3(x3,y第以确定一个圆的方程。设已知圆通过三点P1 (1,3) ,P2(1,-7),P3(6,-2票方程为

(x-x0)2+(y-y)2=R2

按下列步骤求圆的半径R和圆心的坐标(x0,y0)

(1)设圆的一般方程为

x2+y+ax+by+c=0

将圆上三个点的坐标数据分别代入上面方程；

(2) 整理后写关于未知数a,b,c的三元线性方程组AX=b

(3) 用MATLAB解方程组的指令A\ b求解；

(4) 通过配方求出圆的半径R和圆心坐标。

解：(1)将圆上三个点的坐标数据分别代入x2y2a x b y c0可得:

2 _2

1232 a 3b c 0

2

127 a 7b c 0

62 2 2 6a 2b c 0

(2)整理上述方程组得：Ax b

其中，A [1 3 1;1 -7 1;6 -2 1] ； b [-10;-50;-40]

(3) (4)

MATLAB：码：

A=[1 3 1;1 -7 1;6 -2 1];

b=[-10;-50;-40];

x=A\b

r=sqrt((x(1)A2+x(2)A2)/4-x(3)) x0=-(x(1)/2);

y0=-(x(2)/2);

O=[x0,y0]

结果：半径为5,圆心坐标(1, -2)运行结果为: x =

-2

4

-20

r =

5

O =

1 -2