【精品】大连理工矩阵上机作业

【关键字】精品

第一题

Lagrange插值函数

function y=lagrange(x0,y0,x);

n=length(x0);m=length(x);

for i=1:m

z=x(i);

s=0.0;

for k=1:n

p=1.0;

for j=1:n

if j~=k

p=p*(z-x0(j))/(x0(k)-x0(j));

end

end

s=p*y0(k)+s;

end

y(i)=s;

end

x0=[1:10];

y0=[67.052,68.008,69.803,72.024,73.400,72.063,74.669,74.487,74.065,76.777];

lagrange(x0,y0,17)

ans= -1.9516e+12

x=[1:0.1:10];

x=x';

plot(x0,y0,'r');

hold on

plot(x,y,'k');

legend('原函数','拟合函数')

拟合图像如下

拟合函数出现了龙格现象，运用多项式进行插值拟合时，效果并不好，高次多项式会因为误差的不断积累，导致龙格现象的发生。

第二题

function fun =nihe(n)

m=[67.052*10^6,68.008*10^6,69.803*10^6,72.024*10^6,73.400*10^6,72.063*10^6,74.669*10 ^6,74.487*10^6,74.065*10^6,76.777*10^6];

w=[1,2,3,4,5,6,7,8,9,10];

d1=0;d2=0;d3=0;

y1=polyfit(m,w,1);

y2=polyfit(m,w,2);

y3=polyfit(m,w,3);

y2=poly2sym(s2);y3=poly2sym(s3);y4=poly2sym(s4);

f1=subs(y1,17);

f2=subs(y2,17);

f3=subs(y3,17);

for p=1:10;

d1=d1+(subs(y1,w(p))-m(p))^2;

d2=d2+(subs(y2,w(p))-m(p))^2;

d3=d3+(subs(y3,w(p))-m(p))^2;

end

d1=sqrt(d1);

d2=sqrt(d2);

d3=sqrt(d3);

fun=[f1 f2 f3;d2 d3 d4];

return;

结果

三次函数的均方误差最小，拟合的最好。

函数

function f=fun(x)

syms a

a=x;

f=a*a*a+a*a+a-3;

梯度函数

function df=dfun(x)

df=3*x*x+2*x+1;

Newton法

function result=didainewton(x0)

k=0;xk=x0;xi=1;e0=abs(x0-xi);

ek=e0;m=zeros(7,1);

n=zeros(7,1);p=zeros(7,1);

result=zeros(7,3);

while k<7

ak=feval('fun',xk);

bk=feval('dfun',xk);

xk=xk-ak/bk;

e0=ek;

k=k+1;

m(k)=xk;

ek=abs(m(k)-xi);

jingdu=ek/(e0*e0);

n(k)=ek;

p(k)=jingdu;

end

result=[m,n,p];

return;

计算结果

GAUSS消去法

function x=DelGauss(N)

%Gauss??è￥·¨

syms M;

M=N;

a=zeros(M);b=ones(M,1);

for i=1:M

for j=1:M

a(i,j)=1/(i+j-1);

end

end

[n,m]=size(a);nb=length(b); det=1; x=zeros(n,1);

for k=1:n-1

for i=k+1:n

if a(k,k)==0

return

end

m=a(i,k)/a(k,k);

for j=k+1:n

a(i,j)=a(i,j)-m*a(k,j);

end

b(i)=b(i)-m*b(k);

end

det=det*a(k,k);

end

det=det*a(n,n);

for k=n:-1:1

for j=k+1:n

b(k)=b(k)-a(k,j)*x(j);

end

x(k)=b(k)/a(k,k)

end

x=x'

N=4

X=[-4.0000 60.0000 -180.0000 140.0000]T

N=8

X=[ -8.0000001, 504.00001, -7560.0, 46200.0, -138600.0, 216216.0, -168168.0, 51480.0]T

Jacobbi

function x=jacobi(N)

syms M

M=N;

A=zeros(M);b=ones(M,1);

for i=1:M

for j=1:M

A(i,j)=1/(i+j-1);

end

end

p=zeros(M,1);

q=zeros(M,1);

for i=1:M

q(i)=A(i,i);

end

D=diag(q);

L=-tril(A,-1);

U=-triu(A,1);

B=(inv(D))*(L+U);

f=(inv(D))*b;

x=B*p+f;

n=1;

sigmal=1e-8;

while norm((p-x))>=sigmal

p=x;

x=B*p+f;

n=n+1;

end

eval('x');

N=4时

R（B）=2.5821>1，迭代法不收敛

N=8时

R（B）= 6.0422>1,迭代法不收敛

Gauss-Seidel

function x=gaussseidel(N)

syms M

M=N;

A=zeros(M);b=ones(M,1);

for i=1:M

for j=1:M

A(i,j)=1/(i+j-1);

end

end

p=zeros(M,1);

L=-tril(A,-1);

D=diag(diag(A));

U=-triu(A,1);

G=(inv(D-L))*U;

f=(inv(D-L))*b;

x=G*p+f;

n=1;

detal=1e-4;

while abs(x-p)>=detal

p=x;

x=G*p+f;

n=n+1;

end

eval('x');

结果

N=4

X=[2.8527 -14.7133 -3.5388 26.7350]T

N=8

X=[-3.9202 28.0256 -9.4820 -41.3197 -34.0806 -6.4475 28.8462 65.3427]T

共轭梯度法

function result=cg(N)

syms M

M=N;

x0=zeros(M,1);

A=zeros(M);b=ones(M,1);

for i=1:M

for j=1:M

A(i,j)=1/(i+j-1);

end

end

k=0;

r0=b-A*x0;

p0=r0;

sigmal=1e-4;rk=r0;pk=p0;

while dot(r0,r0)>=sigmal

ak=(dot(r0,r0))/(dot(pk,A*pk));

ri=r0;

xk=x0+ak*pk;

rk=r0-ak*A*pk;

bk=(dot(rk,rk))/(dot(ri,ri));

pk=rk+bk*p0;

p0=pk;r0=rk;x0=xk;

k=k+1;

end

result=x0;

return;

结果

N=4

X=[-4.0000 60.0000 -180.0000 140.0000]T

N=8

X=[ 4.5773 -72.4211 199.8408 -9.9954 -174.8581 -171.8511 -10.8652 269.4734]T

第五题

function f=san(x1,y)

n=length(x1);

h=zeros(1,n-1);

for p=1:n-1

h(p)=x1(p+1)-x1(p);

end

for p=1:n-2

r(p)=h(p+1)/(h(p+1)+h(p));

u(p+1)=h(p)/(h(p+1)+h(p));

end

u(1)=1;

r(n-1)=1;

c=2*ones(1,n);

A=diag(r,-1)+diag(u,1)+diag(c);

g=zeros(n,1);

g(1)=3*(y(2)-y(1))/h(1);

g(n)=3*(y(n)-y(n-1))/h(n-1);

for p=1;n-2;

g(p+1)=3*((y(p+2)-y(p+1))/h(p+1)*u(p+1)+r(p)*(y(p+1)-y(p))/h(p)); end

m=A\g;

m=m';

f=m;

return

x1=[0:4];

y=[1,3,3,4,2];

san(x1,y)

ans =

2.5000 1.0000 -0.5000 1.0000 -

3.5000

利用教科书p179的公式（5-35）

function fun=sanci(m,h,xk,y)

syms x

n=length(m);

l=length(xk);

for i=1:n-1

fun=((h(i)+2*(x-xk(i)))*y(i)/h(i)+(x-xk(i))*m(i))*(x-

xk(i+1))^2/h(i)^2+((h(i)-2*(x-xk(i+1)))*y(i+1)/h(i)+(x-

xk(i+1))*m(i+1))*(x-xk(i))^2/h(i)^2

end

end

0<=x<1

((9*x)/2 + 1)*(x - 1)^2 - x^2*(5*x - 8)

1<=x<2

(7*x - 4)*(x - 2)^2 - ((13*x)/2 - 16)*(x - 1)^2

2<=x<3

((11*x)/2 - 8)*(x - 3)^2 - (7*x - 25)*(x - 2)^2

3<=x<4

(9*x - 23)*(x - 4)^2 - ((15*x)/2 - 32)*(x - 3)^2

第六题

函数

function f=fun(x)

syms a;

a=x;

f=a*a*a+2*a*a+10*a-100;

end

梯度函数

function grad=gfun(x)

syms a

a=x;

grad=3*a*a+4*a+10;

迭代函数

f unction result=xianjiefa(x0)

sigmal=1e-6;

maxi=500;

m=zeros(500,1);

d0=fun(x0)/gfun(x0);

x1=x0-d0;

k=1;xk=x1;m(1)=x1;

while k

dk=fun(xk)*(xk-x0)/(fun(xk)-fun(x0)); x0=xk;

xk=xk-dk;

k=k+1;

m(k)=xk;

if abs(xk-x0)

break;

end

end

n=zeros(k,1);

for i=1:k

n(i)=m(i);

end

result=(n);

第七题

function fun=intgra(n)

syms x;

syms m;

sum1=0;

sum2=0;

m=n;

f=(exp(3*x))*cos(pi*x);

x0=0:2*pi/m:2*pi;

a=0;

b=2*pi;

h=(b-a)/m;

for i=1:m-1

xk=a+i*h;

sum1=sum1+subs(f,xk);

end

for i=0:n-1

xk=a+(i+1/2)*h;

sum2=sum2+subs(f,xk);

end

t=((b-a)/(2*n))*(subs(f,a)+2*sum1+subs(f,b));

s=((b-a)/(6*n))*(subs(f,a)+2*sum1+4*sum2+subs(f,b)); z=int(f,x,0,2*pi);

T=vpa(t,2);

S=vpa(s,2);

Z=vpa(z,2)

Dt=abs(T-Z);

Ds=abs(S-Z);

fun=[T;S;Dt;Ds];

return;

计算结果

第八题

Gauss-Chebyshev函数

function result=chebyshev(f,n)

x=zeros(1,n); A=zeros(1,n);

m=n-1;

for i=0:m

x(i+1)=cos((2*i+1)*pi/(2*n));

A(i+1)=pi/n;

end

result=double(sum(A(1:n).*subs(f,x(1:n))));

return

Gauss-Legendre函数

function result=lanendre(f,a,x)

t=pi/4+pi/4*x;

result=double(pi/4*sum(a.*subs(f,t)));

return;

计算过程

syms x

f1=x^2

f2=sin(x)/x

x2=[-0.57735,0.57735]

a2=[1,1];x3=[-0.77460,0,0.77460]

a3=[0.55556,0.88889,0.55556]

x5=[-0.90618,-0.53847,0,0.53847,-0.90618]

a5=[0.23693,0.47863,0.56889,0.47863,0.23693]

fun1=[chebyshev(f1,2),chebyshev(f1,3),chebyshev(f1,5)]

fun2=[lanendre(f2,a2,x2),lanendre(f2,a3,x3),lanendre(f2,a5,x5)] real1=double(int(f1/sqrt(1-x^2),-1,1))

real2=double(int(f2,0,pi/2))

结果

第九题

Euler法

function [t,x]=Euler(f,x0,h)

t=0:h:1;

m=length(t);

x=zeros(1,m);

x(1)=x0;

for i=1:m-1

x(i+1)=x(i)+h*subs(f,{'t','x'},{t(i),x(i)});

end

end

改进Euler法

function [t,x]=ImproveEuler(f,x0,h)

t=0:h:1;

m=length(t);

x=zeros(1,m);

x(1)=x0;

for i=1:m-1

k1=subs(f,{'t','x'},{t(i),x(i)});

x1=x(i)+h*k1;

k2=subs(f,{'t','x'},{t(i+1),x1});

x(i+1)=x(i)+h/2*(k1+k2);

end

end

Runge_Kutta法

function [t,x]=runge_kutta(f,x0,h)

t=0:h:1;

m=length(t);

x=zeros(1,m);

x(1)=x0;

for i=1:m-1

k1=h*subs(f,{'t','x'},{t(i),x(i)});

k2=h*subs(f,{'t','x'},{t(i)+h/2,x(i)+k1/2}); k3=h*subs(f,{'t','x'},{t(i)+h/2,x(i)+k2/2}); k4=h*subs(f,{'t','x'},{t(i)+h,x(i)+k3});

x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;

end

end

调用函数进行计算

syms t x

f=x*cos(t);

x0=1;

h=[0.1,0.01,0.001];

for i=1:3

[t,x1]=Euler(f,x0,h(i));

plot(t,x1);

hold on;

end

>>clc

for i=1:3

[t,x2]=Improved_Euler(f,x0,h(i));

plot(t,x2);

hold on;

end

>>clc

for i=1:3

[t,x3]=Runge_Kutta4(f,x0,h(i));

plot(t,x3);

hold on;

end

此文档是由网络收集并进行重新排版整理.word可编辑版本！