上海大学课程实验报告-数值代数

课程实验报告

COURSE PAPER

课程名称：数值代数与计算方法

课程号：08305114

授课教师：

学号：

姓名：

所属：计算机科学与工程

打印时间：2015

评语：

题目一：

算法：

1、对于问题一：

2、对于问题二：

直接编写递归函数程序，算出三个差分方程的10个近似值程序:

1）

%main.m

clc;

clear all;

a=1;

b=-2;

c=-3;

[x1,x2]=roots(a,b,c)

%roots.m

function [x1,x2]=roots(a,b,c)

d=sqrt(b*b-4*a*c);

if d>0

x1=(-2*c)/(b+d);

x2=(-b-d)/(2*a);

elseif d<0

x1=(-b+d)/(2*a);

x2=(-2*c)/(b-d);

end

2）

%solu.m

function [X,R,P,Q]=solu(X0,R0,P0,P1,Q0,Q1) X(1)=X0;

R(1)=R0;

P(1)=P0;

P(2)=P1;

Q(1)=Q0;

Q(2)=Q1;

for i=1:9

X(i+1)=X(i)/2;

X(i)=X(i+1);

end

for i=1:9

R(i+1)=R(i)/2;

R(i)=R(i+1);

end

for i=3:10

P(i)=3/2*P(i-1)-1/2*P(i-2);

P(i-1)=P(i);

P(i-2)=P(i-1);

end

for i=3:10

Q(i)=5/2*Q(i-1)-Q(i-2);

Q(i-1)=Q(i);

Q(i-2)=Q(i-1);

end

%Xn.m

clc;

clear all;

X0=1;

R0=0.994;

P0=1;

P1=0.497;

Q0=1;

Q1=0.497;

[X,R,P,Q]=solu(X0,R0,P0,P1,Q0,Q1);

for i=1:10

x(i)=i;

end

plot(x,X-R,'r*');

hold on

结果：

1）

x1 =

3

x2 =

-1

2）

题目二：

1、

此题我们以第一个式子为例，求出

g(Xn+1)= (2/(2+3*Xn-Xn.^3))^(1/2)

2、利用二分法：

3、利用牛顿迭代法，并根据题目提示的立方根算法：

此题我们以第三个式子为例：，设A=71/3 ,求出x的3次方便是近似值。

1、

%fixpt.m

function [k,p,err,P]=fixpt(p0,tol,max1)

P(1)=p0;

for k=2:max1

P(k)=g(P(k-1));

err=abs(P(k)-P(k-1));

relerr=err/(abs(P(k))+eps);

p=P(k);

if (err

break;

end

end

if k==max1

disp('maximum number of iterations exceeded')

end

P=P';

%g.m

function y=g(x)

y=sqrt(2/(2+3*x-x.^3));

%iteration.m

clc;

clear all;

p0=0.5;

tol=1.0e-9;

max1=20;

[k,p,err,P]=fixpt(p0,tol,max1);

2、

二分法：

%bisect.m

function [c,err,yc,k,an,cn,bn,ycn]=bisect(a,b,delta) ya=f(a);

yb=f(b);

an(1)=a;

bn(1)=b;

if ya*yb>0

disp('No Root in this district!');

end

max1=1+round((log(b-a)-log(delta))/log(2));

for k=1:max1

c=(a+b)/2;

cn(k)=c;

yc=f(c);

ycn(k)=yc;

if yc==0

a=c;

an(k+1)=a;

b=c;

bn(k+1)=b;

elseif yb*yc>0

b=c;

bn(k+1)=b;

yb=yc;

else

a=c;

an(k+1)=a;

ya=yc;

end

if b-a

break;

end

end

c=(a+b)/2;

cn(max1)=c;

err=abs(b-a);

yc=f(c);

ycn(max1)=yc;

for i=1:9

k(i)=i-1;

end

%f.m

function h=f(x)

h=x*sin(x)-1;

%main.m

clc;

clear all;

a=0;

b=2;

delta=0.000001;

[c,err,yc,k,an,cn,bn,ycn]=bisect(a,b,delta);