数值分析matlab代码

1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6)

disp('牛顿法');

i=1;

n0=180;

p0=pi/3;

tol=10^(-6);

for i=1:n0

p=p0-(p0-sin(p0))/(1-cos(p0));

if abs(p-p0)<=10^(-6)

disp('用牛顿法求得方程的根为')

disp(p);

disp('迭代次数为：')

disp(i)

break;

end

p0=p;

end

if i==n0&&~(abs(p-p0)<=10^(-6))

disp(n0)

disp('次牛顿迭代后无法求出方程的解')

end

2、disp('Steffensen加速');

p0=pi/3;

for i=1:n0

p1=0.5*p0+0.5*cos(p0);

p2=0.5*p1+0.5*cos(p1);

p=p0-((p1-p0).^2)./(p2-2.*p1+p0);

if abs(p-p0)<=10^(-6)

disp('用Steffensen加速求得方程的根为')

disp(p);

disp('迭代次数为：')

disp(i)

break;

end

p0=p;

end

if i==n0&&~(abs(p-p0)<=10^(-6))

disp(n0)

disp('次Steffensen加速后无法求出方程的解')

end

1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，

%误差限为 e=10^-4

disp('二分法')

a=0.2;b=0.26;

tol=0.0001;

n0=10;

fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1;

for i=1:n0

p=(a+b)/2;

fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;

if fp==0||(abs((b-a)/2)

disp('用二分法求得方程的根p=')

disp(p)

disp('二分迭代次数为：')

disp(i)

break;

end

if fa*fp>0

a=p;

else b=p;

end

if i==n0&&~(fp==0||(abs((b-a)/2)

disp(n0)

disp('次二分迭代后没有求出方程的根')

end

2、%使用牛顿法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，

%误差限为 e=10^-4

disp('牛顿法')

p0=0.3;

for i=1:n0

p=p0-(600*(p0.^4)-550*(p0.^3)+200*(p0.^2)-20*p0-1)./(2400*(p0.^3) -1650*p0.^2+400*p0-20);

if(abs(p-p0)

disp('用牛顿法求得方程的根p=')

disp(p)

disp('牛顿迭代次数为：')

disp(i)

break;

end

p0=p;

end

if i==n0&&~(abs(p-p0)

disp(n0)

disp('次牛顿迭代后没有求出方程的根')

end

3、%使用割线法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，

%误差限为 e=10^-4

disp('割线法')

p0=0.2;

p1=0.25;

q0=600*(p0.^4)-550*(p0.^3)+200*(p0.^2)-20*p0-1;

q1=600*(p1.^4)-550*(p1.^3)+200*(p1.^2)-20*p1-1;

for i=2:n0

p=p1-q1*(p1-p0)/(q1-q0);

if abs(p-p1)

disp('用割线法求得方程的根p=')

disp(p)

disp('割线法迭代次数为：')

disp(i)

break;

end

p0=p1;

q0=q1;

pp=p1;

p1=p;

q1=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;

end

if i==n0&&~(abs(p-pp)

disp(n0)

disp('次割线法迭代后没有求出方程的根')

end

4、%使用试位法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，

%误差限为 e=10^-4

disp('试位法')

p0=0.2;

p1=0.25;

q0=600*(p0.^4)-550*(p0.^3)+200*(p0.^2)-20*p0-1;

q1=600*(p1.^4)-550*(p1.^3)+200*(p1.^2)-20*p1-1;

for i=2:n0

p=p1-q1*(p1-p0)/(q1-q0);

if abs(p-p1)

disp('用试位法求得方程的根p=')

disp(p)

disp('试位法迭代次数为：')

disp(i)

break;

end

q=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;

if q*q1<0

p0=p1;

q0=q1;

end

pp=p1;

p1=p;

q1=q;

end

if i==n0&&~(abs(p-pp)

disp(n0)

disp('次试位法迭代后没有求出方程的根')

end

5、%使用muller方法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，

%误差限为 e=10^-4

disp('muller法')

x0=0.1;

x1=0.2;

x2=0.25;

h1=x1-x0;

h2=x2-x1;

d1=((600*(x1.^4)-550*(x1.^3)+200*(x1.^2)-20*x1-1)-(600*(x0.^4)-55 0*(x0.^3)+200*(x0.^2)-20*x0-1))/h1;

d2=((600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)-(600*(x1.^4)-55 0*(x1.^3)+200*(x1.^2)-20*x1-1))/h2;

d=(d2-d1)/(h2+h1);

for i=3:n0

b=d2+h2*d;

D=(b*b-4*(600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)*d)^0.5;

if(abs(d-D)

E=b+D;

else E=b-D;

end

h=-2*(600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)/E;

p=x2+h;

if abs(h)

disp('用muller方法求得方程的根p=')

disp(p)

disp('muller方法迭代次数为：')

disp(i)

break;

end

x0=x1;

x1=x2;

x2=p;

h1=x1-x0;

h2=x2-x1;

d1=((600*(x1.^4)-550*(x1.^3)+200*(x1.^2)-20*x1-1)-(600*(x0.^4)-55 0*(x0.^3)+200*(x0.^2)-20*x0-1))/h1;

d2=((600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)-(600*(x1.^4)-55 0*(x1.^3)+200*(x1.^2)-20*x1-1))/h2;

d=(d2-d1)/(h2+h1);

end

if i==n0%条件有待商榷？！

disp(n0)

disp('次muller方法迭代后没有求出方程的根')

end

1、%观察Lagrange插值的Runge现象

x=-1:0.05:1;

y=1./(1+25.*x.*x);

plot(x,y),grid on;

n=5;

x=-1:2/n:1;

y=1./(1+25.*x.*x);

for i=1:n+1

q(1,i)=y(i);

end

h=0.05;

z=-1:h:1;

for k=1:2/h+1

for i=2:n+1

for j=2:i

q(j,i)=((z(k)-x(i-j+1))*q(j-1,i)-(z(k)-x(i))*q(j-1,i-1))/(x(i)-x( i-j+1));

end

w(k)=q(n+1,n+1);

end

hold on, plot(z,w,'r'),grid on;

%**** n=10 ****

n=10;

x=-1:2/n:1;

y=1./(1+25.*x.*x);

for i=1:n+1

q(1,i)=y(i);

end

h=0.05;

z=-1:h:1;

for k=1:2/h+1

for i=2:n+1

for j=2:i

q(j,i)=((z(k)-x(i-j+1))*q(j-1,i)-(z(k)-x(i))*q(j-1,i-1))/(x(i)-x( i-j+1));

end

w(k)=q(n+1,n+1);

end

hold on,plot(z,w,'k'),grid on;

legend ('原始图','n=5','n=10');

2、%固支样条插植

%********第一段********

x=[1,2,5,6,7,8,10,13,17];

a=[3,3.7,3.9,4.2,5.7,6.6,7.1,6.7,4.5];

n=numel(a);

for i=1:n-1

h(i)=x(i+1)-x(i);

end

A=[2*h(1),h(1),0,0,0,0,0,0,0;

h(1),2*(h(1)+h(2)),h(2),0,0,0,0,0,0;

0,h(2),2*(h(2)+h(3)),h(3),0,0,0,0,0;

0,0,h(3),2*(h(3)+h(4)),h(4),0,0,0,0;

0,0,0,h(4),2*(h(4)+h(5)),h(5),0,0,0;

0,0,0,0,h(5),2*(h(5)+h(6)),h(6),0,0;

0,0,0,0,0,h(6),2*(h(6)+h(7)),h(7),0;

0,0,0,0,0,0,h(7),2*(h(7)+h(8)),h(8);

0,0,0,0,0,0,0,h(8),2*h(8)];

e=[3*(a(2)-a(1))/h(1)-3;

3*(a(3)-a(2))/h(2)-3*(a(2)-a(1))/h(1);

3*(a(4)-a(3))/h(3)-3*(a(3)-a(2))/h(2);

3*(a(5)-a(4))/h(4)-3*(a(4)-a(3))/h(3);

3*(a(6)-a(5))/h(5)-3*(a(5)-a(4))/h(4);

3*(a(7)-a(6))/h(6)-3*(a(6)-a(5))/h(5);

3*(a(8)-a(7))/h(7)-3*(a(7)-a(6))/h(6);

3*(a(9)-a(8))/h(8)-3*(a(8)-a(7))/h(7);

3*(-0.67)-3*(a(9)-a(8))/h(8)];

c=inv(A)*e;

for i=1:8

b(i)=(a(i+1)-a(i))/h(i)-h(i)*(2*c(i)+c(i+1))/3;

d(i)=(c(i+1)-c(i))/(3*h(i));

end

for i=1:8

z=x(i):0.05:x(i+1);

w=a(i)+b(i).*(z-x(i))+c(i).*(z-x(i)).^2+d(i).*(z-x(i)).^3; grid on, plot(z,w),hold on;

end

%********第二段********

x=[17,20,23,24,25,27,27.7];

a=[4.5,7,6.1,5.6,5.8,5.2,4.1];

for i=1:6

h(i)=x(i+1)-x(i);

end

A=[2*h(1),h(1),0,0,0,0,0;

h(1),2*(h(1)+h(2)),h(2),0,0,0,0;

0,h(2),2*(h(2)+h(3)),h(3),0,0,0;

0,0,h(3),2*(h(3)+h(4)),h(4),0,0;

0,0,0,h(4),2*(h(4)+h(5)),h(5),0;

0,0,0,0,h(5),2*(h(5)+h(6)),h(6)

0,0,0,0,0,h(6),2*h(6)];

e=[3*(a(2)-a(1))/h(1)-3*3;

3*(a(3)-a(2))/h(2)-3*(a(2)-a(1))/h(1);

3*(a(4)-a(3))/h(3)-3*(a(3)-a(2))/h(2);

3*(a(5)-a(4))/h(4)-3*(a(4)-a(3))/h(3);

3*(a(6)-a(5))/h(5)-3*(a(5)-a(4))/h(4);

3*(a(7)-a(6))/h(6)-3*(a(6)-a(5))/h(5);

3*(-4)-3*(a(7)-a(6))/h(6)];

c=inv(A)*e;

for i=1:6

b(i)=(a(i+1)-a(i))/h(i)-h(i)*(2*c(i)+c(i+1))/3;

d(i)=(c(i+1)-c(i))/(3*h(i));

end

for i=1:6

z=x(i):0.05:x(i+1);

w=a(i)+b(i).*(z-x(i))+c(i).*(z-x(i)).^2+d(i).*(z-x(i)).^3; grid on, plot(z,w),hold on;

end

%********第三段********

x=[27.7,28,29,30];

a=[4.1,4.3,4.1,3];

for i=1:3

h(i)=x(i+1)-x(i);

end

A=[2*h(1),h(1),0,0;

h(1),2*(h(1)+h(2)),h(2),0;

0,h(2),2*(h(2)+h(3)),h(3);

0,0,h(3),2*h(3)];

e=[3*(a(2)-a(1))/h(1)-3*0.33;

3*(a(3)-a(2))/h(2)-3*(a(2)-a(1))/h(1);

3*(a(4)-a(3))/h(3)-3*(a(3)-a(2))/h(2);

3*(-1.5)-3*(a(4)-a(3))/h(3)];

c=inv(A)*e;

for i=1:3

b(i)=(a(i+1)-a(i))/h(i)-h(i)*(2*c(i)+c(i+1))/3;

d(i)=(c(i+1)-c(i))/(3*h(i));

end

for i=1:3

z=x(i):0.05:x(i+1);

w=a(i)+b(i).*(z-x(i))+c(i).*(z-x(i)).^2+d(i).*(z-x(i)).^3; grid on, plot(z,w),hold on;

end

grid on,title('注：横纵坐标的比例不一样！！！');

1、%用不动点迭代法求方程 x-e^x+4=0的正根与负根，误差限是10^-6%

disp('不动点迭代法');

n0=100;

p0=-5;

for i=1:n0

p=exp(p0)-4;

if abs(p-p0)<=10^(-6)

if p<0

disp('|p-p0|=')

disp(abs(p-p0))

disp('不动点迭代法求得方程的负根为:')

disp(p);

break;

else

disp('不动点迭代法无法求出方程的负根.')

end

else

p0=p;

end

if i==n0

disp(n0)

disp('次不动点迭代后无法求出方程的负根')

end

p1=1.7;

for i=1:n0

pp=exp(p1)-4;

if abs(pp-p1)<=10^(-6)

if pp>0

disp('|p-p1|=')

disp(abs(pp-p1))

disp('用不动点迭代法求得方程的正根为')

disp(pp);

else

disp('用不动点迭代法无法求出方程的正根');

end

break;

else

p1=pp;

end

if i==n0

disp(n0)

disp('次不动点迭代后无法求出方程的正根')

end

2、%用牛顿法求方程 x-e^x+4=0的正根与负根，误差限是10^-6 disp('牛顿法')

n0=80;

p0=1;

for i=1:n0

p=p0-(p0-exp(p0)+4)/(1-exp(p0));

if abs(p-p0)<=10^(-6)

disp('|p-p0|=')

disp(abs(p-p0))

disp('用牛顿法求得方程的正根为')

disp(p);

break;

else

p0=p;

end

if i==n0

disp(n0)

disp('次牛顿迭代后无法求出方程的解')

end

p1=-3;

for i=1:n0

p=p1-(p1-exp(p1)+4)/(1-exp(p1));

if abs(p-p1)<=10^(-6)

disp('|p-p1|=')

disp(abs(p-p1))

disp('用牛顿法求得方程的负根为')

disp(p);

break;

else

p1=p;

end

if i==n0

disp(n0)

disp('次牛顿迭代后无法求出方程的解')

end

1、使用欧拉法、改进欧拉法和四阶R-K方法求下列微分方程的解。

≤

)(

0,1

'。其精确解为

)0(

y t+

≤

%欧拉法

disp('欧拉法');

a=0;

b=5;

A=1;

h=0.05;

for k=1:3

t=a;

w=A;

for i=1:(b-a)/h

w=w+h*(-w+t+1);

t=a+i*h;

end

disp('h=');

disp(h);

disp('计算结果是');

disp(w);

disp('误差为');

disp(abs(exp(-5)+5-w));

h=h*2;

end

%改进欧拉法

disp('改进欧拉法');

a=0;

b=5;

A=1;

h=0.05;

for k=1:3

t=a;

w=A;

for i=1:(b-a)/h

w=w+h*((-w+t+1)+(-((-w+t+1)*h+w)+(t+h)+1))/2; t=a+i*h;

end

disp('h=');

disp(h);

disp('计算结果是');

disp(w);

disp('误差为');

disp(abs(exp(-5)+5-w));

h=h*2;

end

%四阶R-K方法

disp('四阶R-K方法');

a=0;

b=5;

A=1;

h=0.05;

for k=1:3

t=a;

w=A;

for i=1:(b-a)/h

k1=h*(-w+t+1);

k2=h*(-(w+k1/2)+(t+h/2)+1);

k3=h*(-(w+k2/2)+(t+h/2)+1);

k4=h*(-(w+k3)+(t+h)+1);

w=w+(k1+2*k2+2*k3+k4)/6;

t=a+i*h;

end

disp('h=');

disp(h);

disp('计算结果是');

disp(w);

disp('误差为');

disp(abs(exp(-5)+5-w));

h=h*2;

end

1、题目：使用Romberg 积分，对于

?+4802)cos (1dx x 计算下列各近似值。

a.确定1,51,41,31,21,1,,,,R R R R R

b.确定5,54,43,32,2,,,R R R R

c.6,65,64,63,62,61,6,,,,,R R R R R R

d.确定10,109,98,87,7,,,R R R R

n=5;

a=0;

b=48;

h(1,1)=b-a;

disp('-------求R11，R21，R31，R41，R51-----------');

fa=sqrt(1+(cos(a))^2);

fb=sqrt(1+(cos(b))^2);

r(1,1)=h(1,1)/2*(fa+fb);

disp('R11，R21，R31，R41，R51分别为');

disp(r(1,1));

for i=2:n

h(i,1)=(b-a)/(2^(i-1));

sum=0;

for k=1:2^(i-2)

x=a+(2*k-1)*h(i,1);

sum=sum+sqrt(1+(cos(x)).^2);

end

r(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);

disp(r(i,1));

end

disp('-----求R22,R33，R44,R55---------');

disp('R22，R33，R44，R55分别为');

for k=2:n

for j=2:k

r(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);

end

disp(r(k,k));

end

disp('-----R61，R62，R63，R64，R65，R66-----');

disp('R61，R62，R63，R64，R65，R66分别为');

n=6;

for i=2:n

h(i,1)=(b-a)/(2^(i-1));

sum=0;

for k=1:2^(i-2)

x=a+(2*k-1)*h(i,1);

sum=sum+sqrt(1+(cos(x)).^2);

end

r(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);

end

for k=2:n

for j=2:k

r(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);

end

for i=1:n

disp(r(6,i));

end

disp('-----R77，R88，R99，R10,10--------');

disp('R77，R88，R99，R10,10分别为');

n=10;

for i=2:n

h(i,1)=(b-a)/(2^(i-1));

sum=0;

for k=1:2^(i-2)

x=a+(2*k-1)*h(i,1);

sum=sum+sqrt(1+(cos(x)).^2);

end

r(i,1)=0.5*(r(i-1,1)+h(i-1,1)*sum);

end

for k=2:n

for j=2:k

r(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);

end

for i=7:10

disp(r(i,i));

end

disp('为什么龙贝格积分的结果是不同的？');

disp('因为有舍入误差和截断误差的影响');

1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6)

disp('牛顿法')

i=1;

n0=180;

p0=pi/3;

tol=10^(-6);

for i=1:n0

p=p0-(p0-sin(p0))/(1-cos(p0));

if abs(p-p0)<=10^(-6)

%disp('|p-p0|=')

%disp(abs(p-p0))

disp('用牛顿法求得方程的根为')

disp(p);

disp('迭代次数为：')

disp(i)

break;

end

p0=p;

end

if i==n0

disp(n0)

disp('次计算后无法求出方程的解')

end

2、%用求重根的方法求f(x)的零点

disp('求重根方法一');

p0=pi/3;

for i=1:n0 p=p0-((p0-sin(p0)).*(1-cos(p0)))./((1-cos(p0)).^2-(p0-sin(p0)).*s

in(p0));

if(abs(p-p0)

disp('用求重根的方法一求得方程的根p=')

disp(p)

disp('迭代次数为：')

disp(i)

break;

else

p0=p;

end

if i==n0&&~(abs(p-p0)

disp(n0)

disp('次迭代后没有用求重根的方法求出方程的根')

end

3、%用求重根的方法求f(x)的零点

disp('求重根方法二');

p0=pi/3;

for i=1:n0

p=p0-1*(p0-sin(p0))/(1-cos(p0));

if (abs(p-p0)

disp('用求重根的方法二求得方程的根p=')

disp(p)

disp('迭代次数为：')

disp(i)

break ;

end

p0=p;

end

if i==n0&&~(abs(p-p0)

disp(n0)

disp('次迭代后没有求出方程的根')

end

1、题目：b

Ax Crout ,9)2.i (01b 10b 1-,2a ,9)2,i (1-,210101019,101210,10111-,1,=??========??====?+分解，求解使用三对角方程组给出。对每个和维列向量，由是给出。设和对每个三对角矩阵，由是设i i i i i ii b b a a a a a a A

disp('用三对角方程组Crout 分解，求解方程组');

for i=1:10

v1(i)=2;

end

for i=1:9

v2(i)=-1;

v3(i)=-1;

end

A=diag(v1,0)+diag(v2,1)+diag(v3,-1);

c(1)=1;

c(10)=1;

for i=2:9

c(i)=0;

end

for i=1:10

A(i,11)=c(i);

end

n=10;

l(1,1)=A(1,1);

u(1,2)=A(1,2)/l(1,1);

z(1)=A(1,n+1)/l(1,1);

for i=2:(n-1)

l(i,i-1)=A(i,i-1);

l(i,i)=A(i,i)-l(i,i-1)*u(i-1,i);

u(i,i+1)=A(i,i+1)/l(i,i);

z(i)=(A(i,n+1)-l(i,i-1)*z(i-1))/l(i,i);

end

l(n,n-1)=A(n,n-1);

l(n,n)=A(n,n)-l(n,n-1)*u(n-1,n);

z(n)=(A(n,n+1)-l(n,n-1)*z(n-1))/l(n,n);

x(n)=z(n);

for i=(n-1):(-1):1

x(i)=z(i)-u(i,i+1)*x(i+1);

end

disp('方程组的解X=');

disp(x);

题目：用四阶R-K 方法求下列初值问题的解。

1、

t t t t t t t e t e e t u e e e t u h u t e t t u u u u t e t u u u 22-522-512222*********

231)(31-31)(,2.01

)0(.10,)4-2(-4'1

)0(.10,)12(-23'++=+===≤≤++==≤≤++=和精确解为

2、 2--2-6

1)(1

.0,0)0('y(0),1,0-'2'-'3t e te e t t y h y t t te y y y t t t t +====≤≤=+精确解为 1、%用四阶R-K 方法求P322的1a

disp('P322 1(a)');

a=0;

b=1;

h=0.02;

n=(b-a)/h;

t=a;

u1=1;

u2=1;

for i=1:n

%每个方程的k1

k(1,1)=h*(3*u1+2*u2-(2*t^2+1)*exp(2*t));

k(1,2)=h*(4*u1+u2+(t^2+2*t-4)*exp(2*t));

%每个方程的k2

k(2,1)=h*(3*(u1+0.5*k(1,1))+2*(u2+0.5*k(1,2))-(2*(t+0.5*h)^2+1)*e xp(2*(t+0.5*h)));

k(2,2)=h*(4*(u1+0.5*k(1,1))+(u2+0.5*k(1,2))+((t+0.5*h)^2+2*(t+0.5 *h)-4)*exp(2*(t+0.5*h)));

%每个方程的k3

k(3,1)=h*(3*(u1+0.5*k(2,1))+2*(u2+0.5*k(2,2))-(2*(t+0.5*h)^2+1)*e xp(2*(t+0.5*h)));

k(3,2)=h*(4*(u1+0.5*k(2,1))+(u2+0.5*k(2,2))+((t+0.5*h)^2+2*(t+0.5 *h)-4)*exp(2*(t+0.5*h)));

%每个方程的k4

k(4,1)=h*(3*(u1+k(3,1))+2*(u2+k(3,2))-(2*(t+h)^2+1)*exp(2*(t+h))) ;

k(4,2)=h*(4*(u1+k(3,1))+(u2+k(3,2))+((t+h)^2+2*(t+h)-4)*exp(2*(t+ h)));

%求每个方程的w

u1=u1+(k(1,1)+2*k(2,1)+2*k(3,1)+k(4,1))/6;

u2=u2+(k(1,2)+2*k(2,2)+2*k(3,2)+k(4,2))/6;

t=a+i*h;

end

disp('w1=');

disp(u1);

disp('误差为');

disp(abs(u1-(exp(5)/3-exp(-1)/3+exp(2))));

disp('w2=');

disp(u2);

disp('误差为');

disp(abs(u2-exp(5)/3-exp(-1)/3*2-exp(2)));

2、%用四阶R-K方法求P322的2a

disp('P322 2(a)');

a=0;

b=1;

h=0.1;

n=(b-a)/h;

t=a;

u1=0;

u2=0;

for i=1:n

%每个方程的k1

k(1,1)=h*(u2);

k(1,2)=h*(-u1+2*u2+t*exp(t)-t);

%每个方程的k2

k(2,1)=h*(u2+0.5*k(1,2));

k(2,2)=h*(-(u1+0.5*k(1,1))+2*(u2+0.5*k(1,2))+(t+0.5*h)*exp(t+0.5*h)-(t+0.5*h));

%每个方程的k3

k(3,1)=h*(u2+0.5*k(2,2));

k(3,2)=h*(-(u1+0.5*k(2,1))+2*(u2+0.5*k(2,2))+(t+0.5*h)*exp(t+0.5*h)-(t+0.5*h));

%每个方程的k4

k(4,1)=h*(u2+k(3,2));

k(4,2)=h*(-(u1+k(3,1))+2*(u2+k(3,2))+(t+h)*exp(t+h)-(t+h)); %求每个方程的w

u1=u1+(k(1,1)+2*k(2,1)+2*k(3,1)+k(4,1))/6;

u2=u2+(k(1,2)+2*k(2,2)+2*k(3,2)+k(4,2))/6;

t=a+i*h;

end

disp('y1=');

disp(u1);

disp('误差为');

disp(abs(u1-exp(1)/6-exp(1)+3));

题目：分别用高斯消元法、选主元高斯消元法和部分选主元高斯消元法解下列方程组：

)

18260870.1,-02065405.0,-01269269.0,17682530.0(16

.4-1.13-110.03.1515

.29.99-10044

.3-2.12122.0-2.1412

.1100-11.219.1432143243214321精确解为=+=+=+=++x x x x x x x x x x x x x x x

1、disp('高斯消元法');

n=4;

a=[1.19,2.11,-100,1,1.12;

14.2,-0.122,12.2,-1,3.44;

0,100,-99.9,1,2.15;

15.3,0.11,-13.1,-1,4.16];

p=1;

for i=1:(n-1)

for j=i:n

if a(j,i)~=0

p=j;

break;

end

if a(j,i)==0

disp('该方程组没有唯一解');

end

if p~=i

b=a(i,:);

a(i,:)=a(p,:);

a(p,:)=b;

end

for j=(i+1):n

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

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

end

if a(n,n)==0

disp('该方程组没有唯一解');

end

x(n)=a(n,n+1)/a(n,n);

for i=(n-1):(-1):1

sum=0;

for j=(i+1):n

sum=sum+a(i,j)*x(j);

end

x(i)=(a(i,n+1)-sum)/a(i,i); end

disp('方程组的解 X=');

disp(x);

2、disp('部分选主元高斯消元法');

a=[1.19,2.11,-100,1,1.12;

14.2,-0.122,12.2,-1,3.44;

0,100,-99.9,1,2.15;

15.3,0.11,-13.1,-1,4.16];

for i=1:(n-1)

b=abs(a(i,i));

for j=i:n

if abs(a(j,i))>b

p=j;

end

if a(p,i)==0

disp('该方程组没有唯一解');

end

if i~=p

b=a(i,:);

a(i,:)=a(p,:);

a(p,:)=b;

end

for j=(i+1):n

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

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

end

if a(n,n)==0

disp('该方程组没有唯一解'); end

x(n)=a(n,n+1)/a(n,n);

for i=(n-1):(-1):1

sum=0;

for j=(i+1):n

sum=sum+a(i,j)*x(j);

end

x(i)=(a(i,n+1)-sum)/a(i,i); end

disp('方程组的解 X=');

disp(x);

3、disp('按比例选主元高斯消元法'); a=[1.19,2.11,-100,1,1.12;

14.2,-0.122,12.2,-1,3.44;

0,100,-99.9,1,2.15;

15.3,0.11,-13.1,-1,4.16];

for i=1:n

s(i)=abs(a(i,1));

for j=1:n

if abs(a(i,j))>s(i)

s(i)=abs(a(i,j));

end

if s(i)==0

disp('该方程组没有唯一解');

break;

end

for i=1:(n-1)

《MATLAB与数值分析》第一次上机实验报告

电子科技大学电子工程学院标准实验报告（实验）课程名称MATLAB与数值分析学生姓名：李培睿学号：2013020904026 指导教师：程建

一、实验名称《MATLAB与数值分析》第一次上机实验二、实验目的 1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算操作。（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序） 2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。（用.m文件编写进行符号因式分解和函数求反的程序） 3. 掌握Matlab函数的编写规范。 4、掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、三维曲线和面的填充、三维等高线等。（用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释） 5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。三、实验内容 1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。并以x， y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图，允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。不允许使用sort 函数。四、实验数据及结果分析题目一： ①在Editor窗口编写函数代码如下：

数值分析MATLAB上机实验

数值分析实习报告姓名：gestepoA 学号：201******* 班级：***班

序言随着计算机技术的迅速发展，数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。要解决工程问题，往往需要处理很多数学模型，不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：解法是否稳定、是否收敛及熟练的速度等。而且还能减少大量的人工计算。由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如MATLAB，C++，VB，JAVA的辅助软件来解决，得到一个满足误差限的解。本文所计算题目，均采用MATLAB进行编程，MATLAB被称为第四代计算机语言，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁，它用更直观的、符合人们思维习惯的代码。它具有以下优点： 1友好的工作平台和编程环境。MATLAB界面精致，人机交互性强，操作简单。 2简单易用的程序语言。MATLAB是一个高级的矩阵/阵列语言，包含控制语言、函数、数据结构，具有输入、输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（M 文件）后再一起运行。 3强大的科学计算机数据处理能力。包含大量计算算法的集合，拥有600多个工程中要用到的数学运算函数。 4出色的图像处理功能，可以方便地输出二维图像，便于我们绘制函数图像。

目录 1 第一题 (4) 1.1 实验目的 (4) 1.2 实验原理和方法 (4) 1.3 实验结果 (5) 1.3.1 最佳平方逼近法 (5) 1.3.2 拉格朗日插值法 (7) 1.3.3 对比 (8) 2 第二题 (9) 2.1实验目的 (9) 2.2 实验原理和方法 (10) 2.3 实验结果 (10) 2.3.1 第一问 (10) 2.3.2 第二问 (11) 2.3.3 第三问 (11) 3 第三题 (12) 3.1实验目的 (12) 3.2 实验原理和方法 (12) 3.3 实验结果 (12) 4 MATLAB程序 (14)

matlab数值计算(命令与示例)

MATLAB数值计算 MATLAB数值计算 (1) 1创建矩阵 (3) 1.1直接输入 (3) 1.2向量 (3) 1.2.1linspace：线性分布 (3) 1.2.2冒号法 (3) 1.3函数创建 (4) 1.3.1eye:单位矩阵 (4) 1.3.2rand：随机矩阵 (4)

1.3.3zeros：全0矩阵 (4) 1.3.4ones：全1矩阵 (5) 2矩阵运算 (5) 2.1加减 (5) 2.1.1[M×N]±[M×N] (5) 2.2乘 (6) 2.2.1[M×N]*a (6) 2.2.2[M×N]*[N×M] (6) 2.3乘方 (7) 2.3.1[M×M]^a (7) 2.3.2a^[M×M] (7) 2.4特殊运算 (8) 2.4.1求逆inv (8) 2.4.2行列式det (8) 2.4.3特征值eig (8) 2.4.4转置'和.' (9) 2.4.5变形reshape (10) 2.4.6翻转rot90,fliplr,flipud (11) 2.4.7抽取diag,tril,triu (12) 2.5数组运算 (12) 2.5.1乘 (12) [M×N].*[M×N] (12) 2.5.2除 (13) [M×N]./[M×N] (14) [M×N].\[M×N] (14) 2.5.3乘方 (14) [M×N].^[M×N] (15) a.^[M×N] (15) 2.6除法 (15) 2.6.1求解线性方程组 (15) 3多项式 (16) 3.1系数表示法poly (16) 3.2求根roots (16) 3.3乘法conv (16) 3.4除法deconv (17) 3.5求值polyval (17) 3.6微分polyder (18)

Matlab作业3(数值分析)答案

Matlab作业3（数值分析）机电工程学院（院、系）专业班组学号姓名实验日期教师评定 1.计算多项式乘法(x2+2x+2)(x2+5x+4)。答： 2. （1）将(x-6)(x-3)(x-8)展开为系数多项式的形式。（2）求解在x=8时多项式(x-1)(x-2) (x-3)(x-4)的值。答：（1）（2）

3. y=sin(x)，x从0到2π，?x=0.02π，求y的最大值、最小值、均值和标准差。 4.设ｘ＝[0.00.30.8 1.1 1.6 2.3]'，y=[0.500.82 1.14 1.25 1.35 1.40]'，试求二次多项式拟合系数，并据此计算x1=[0.9 1.2]时对应的y1。解：x=[0.0 0.3 0.8 1.1 1.6 2.3]'; %输入变量数据x y=[0.50 0.82 1.14 1.25 1.35 1.40]'; %输入变量数据y p=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数p x1=[0.9 1.2]; %输入点x1 y1=polyval(p,x1) %估计x1处对应的y1 p = -0.2387 0.9191 0.5318 y1 = a) 1.2909

5.实验数据处理：已知某压力传感器的测试数据如下表 p为压力值，u为电压值，试用多项式 d cp bp ap p u+ + + =2 3 ) ( 来拟合其特性函数，求出a,b,c,d，并把拟合曲线和各个测试数据点画在同一幅图上。解： >> p=[0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9]; u=[10,11,13,14,17,18,22,24,29,34,39]; x=polyfit(p,u,3) %得多项式系数 t=linspace(0,10,100); y=polyval(x,t); %求多项式得值 plot(p,u,'*',t,y,'r') %画拟和曲线 x = 0.0195 -0.0412 1.4469 9.8267

数值分析的matlab实现

第2章牛顿插值法实现参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55. 求牛顿插值多项式和差商的MA TLAB 主程序： function[A,C,L,wcgs,Cw]=newpoly(X,Y) n=length(X);A=zeros(n,n);A(:,1) =Y'; s=0.0;p=1.0;q=1.0;c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1; end C=A(n,n);b=poly(X(n));q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C);Q=poly2sym(q1); syms M wcgs=M*Q/c1;Cw=q1/c1; （1）保存名为newpoly.m 的M 文件（2）输入MA TLAB 程序 >> X=[242,243,249,250]; >> Y=[13.681,13.526,13.098,13.095]; >> [A,C,L,wcgs,Cw]=newpoly(X,Y) 输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其 ) ()()1(ξ+n n f x R 的系数向量Cw 。 A = 13.6810 0 0 0 13.5260 -0.1550 0 0 13.0980 -0.0713 0.0120 0 13.0950 -0.0030 0.0098 -0.0003 C = 1.0e+003 *

matlab与数值分析作业

数值分析作业(1) 1：思考题（判断是否正确并阐述理由）（a）一个问题的病态性如何，与求解它的算法有关系。 (b)无论问题是否病态，好的算法都会得到它好的近似解。 (c)计算中使用更高的精度，可以改善问题的病态性。 (d)用一个稳定的算法计算一个良态问题，一定会得到他好的近似解。 (e)浮点数在整个数轴上是均匀分布。 (f)浮点数的加法满足结合律。 (g)浮点数的加法满足交换律。 (h)浮点数构成有效集合。 (i)用一个收敛的算法计算一个良态问题，一定得到它好的近似解。√2: 解释下面Matlab程序的输出结果 t=0.1; n=1:10; e=n/10-n*t 3：对二次代数方程的求解问题 20 ++= ax bx c 有两种等价的一元二次方程求解公式

2224b x a c x b ac -±==- 对 a=1，b=-100000000，c=1，应采用哪种算法？ 4：函数sin x 的幂级数展开为： 357 sin 3!5!7! x x x x x =-+-+ 利用该公式的Matlab 程序为 function y=powersin(x) % powersin. Power series for sin(x) % powersin(x) tries to compute sin(x)from a power series s=0; t=x; n=1; while s+t~=s; s=s+t; t=-x^2/((n+1)*(n+2))*t n=n+2; end

（a ）解释上述程序的终止准则；（b ）对于x=/2π、x=11/2π、x =21/2π，计算的精度是多少？分别需要计算多少项？ 5:指数函数的幂级数展开 2312!3!x x x e x =+++ + 根据该展开式，编写Matlab 程序计算指数函数的值，并分析计算结果（重点分析0x <的计算结果）。

MATLAB与数值分析课程总结

MATLAB与数值分析课程总结姓名：董建伟学号：2015020904027 一：MATLAB部分 1.处理矩阵－容易矩阵的创建（1）直接创建注意 a中括号里可以用空格或者逗号将矩阵元素分开 b矩阵元素可以是任何MATLAB表达式，如实数复数等 c可以调用赋值过的任何变量，变量名不要重复，否则会被覆盖（2）用MATLAB函数创建矩阵如：a空阵［］ b rand/randn——随机矩阵 c eye——单位矩阵 d zeros ——0矩阵 e ones——1矩阵 f magic——产生n阶幻方矩阵等向量的生成（1）用冒号生成向量（2）使用linspace和logspace分别生成线性等分向量和对数等分向量矩阵的标识和引用（1）向量标识（2）“0 1”逻辑向量或矩阵标识（3）全下标，单下标，逻辑矩阵方式引用字符串数组（1）字符串按行向量进行储存（2）所有字符串用单引号括起来（3）直接进行创建矩阵运算（1）注意与数组点乘，除与直接乘除的区别，数组为乘方对应元素的幂

（2）左右除时斜杠底部靠近谁谁是分母（3）其他运算如，inv矩阵求逆，det行列式的值, eig特征值，diag 对角矩阵 2.绘图－轻松 plot－绘制二维曲线（1）plot(x)绘制以x为纵坐标的二维曲线 plot(x,y) 绘制以x为横坐标，y为纵坐标的二维曲线 x,y为向量或矩阵（2）plot（x1,y1,x2,y2，。。。。。。）绘制多条曲线，不同字母代替不同颜色:b蓝色，y黄色，r红色，g绿色（3）hold on后面的pl ot图像叠加在一起 hold off解除hold on命令，plot将先冲去窗口已有图形（4）在hold后面加上figure，可以绘制多幅图形（5）subplot在同一窗口画多个子图三维图形的绘制 (1)plot3(x,y,z,’s’) s是指定线型，色彩，数据点形的字符串 (2)［X,Y］＝meshgrid(x,y)生成平面网格点 (3)mesh(x,y,z，c)生成三维网格点，c为颜色矩阵 (4)三维表面处理mesh命令对网格着色，surf对网格片着色 (5)contour绘制二维等高线 (6)axis(［x1,xu,y1,yu］)定义x,y的显示范围 3．编程－简洁（1）变量命名时可以由字母，数字，下划线，但是不得包含空格和标点（2）最常用的数据类型只有双精度型和字符型，其他数据类型只在特殊条件下使用（3）为得到高效代码，尽量提高代码的向量化程度，避免使用循环结构

数值分析的MATLAB程序

列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1);

end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x')

matlab实现数值分析插值及积分

Matlab实现数值分析插值及积分摘要：数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中，常常将实际问题转化为数学模型来解决，这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要计算差值和积分，对于问题一，可以分别利用朗格朗日插值公式，牛顿插值公式，埃特金逐次线性插值公式来进行编程求解，具体matlab代码见正文。编程求解出来的结果为：=+。其中Aitken插值计算的结果图如下：对于问题二，可以分别利用复化梯形公式，复化的辛卜生公式，复化的柯特斯公式编写程序来进行求解，具体matlab代码见正文。编程求解出来的结果为： 0.6932 其中复化梯形公式计算的结果图如下：

问题重述问题一：已知列表函数表格 1 分别用拉格朗日，牛顿，埃特金插值方法计算。问题二：用复化的梯形公式，复化的辛卜生公式，复化的柯特斯公式计算积分，使精度小于5。问题解决问题一：插值方法对于问题一，用三种差值方法：拉格朗日，牛顿，埃特金差值方法来解决。一、拉格朗日插值法：拉格朗日插值多项式如下：首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数，对任一点i x 所对应的插值基函数 )(x l i ，由于在所有),,1,1,,1,0(n i i j x j +-=取零值，因此)(x l i 有因子 )())(()(110n i i x x x x x x x x ----+- 。又因)(x l i 是一个次数不超过n 的多项式，所以只可能相差一个常数因子，固)(x l i 可表示成： )())(()()(110n i i i x x x x x x x x A x l ----=+- 利用1)(=i i x l 得：

同济大学数值分析matlab编程题汇编

MATLAB 编程题库 1.下面的数据表近似地满足函数2 1cx b ax y ++=，请适当变换成为线性最小二乘问题，编程求最好的系数c b a ,,，并在同一个图上画出所有数据和函数图像. 625 .0718.0801.0823.0802.0687.0606.0356.0995 .0628.0544.0008.0213.0362.0586.0931.0i i y x ---- 解： x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; A=[x ones(8,1) -x.^2.*y]; z=A\y; a=z(1); b=z(2); c=z(3); xh=-1:0.1:1; yh=(a.*xh+b)./(1+c.*xh.^2); plot(x,y,'r+',xh,yh,'b*')

2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数 10 的近似根，并写出调用方式：精度为10 解： >> edit gexianfa.m function [x iter]=gexianfa(f,x0,x1,tol) iter=0; while(norm(x1-x0)>tol) iter=iter+1; x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end >> edit f.m function v=f(x) v=x.*log(x)-1; >> edit g.m function z=g(y) z=y.^5+y-1; >> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 = 1.7632 iter1 = 6 >> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 = 0.7549 iter2 = 8

第2讲 matlab的数值分析

第二讲MATLAB的数值分析 2-1矩阵运算与数组运算矩阵运算和数组运算是MATLAB数值运算的两大类型，矩阵运算是按矩阵的运算规则进行的，而数组运算则是按数组元素逐一进行的。因此，在进行某些运算（如乘、除）时，矩阵运算和数组运算有着较大的差别。在MATLAB中，可以对矩阵进行数组运算，这时是把矩阵视为数组，运算按数组的运算规则。也可以对数组进行矩阵运算，这时是把数组视为矩阵，运算按矩阵的运算规则进行。 1、矩阵加减与数组加减矩阵加减与数组加减运算效果一致，运算符也相同，可分为两种情况：（1）若参与运算的两矩阵（数组）的维数相同，则加减运算的结果是将两矩阵的对应元素进行加减，如 A=[1 1 1;2 2 2;3 3 3]; B=A; A+B ans= 2 2 2 4 4 4 6 6 6 （2）若参与运算的两矩阵之一为标量（1*1的矩阵），则加减运算的结果是将矩阵（数组）的每一元素与该标量逐一相加减，如 A=[1 1 1;2 2 2;3 3 3]; A+2 ans= 3 3 3 4 4 4 5 5 5 2、矩阵乘与数组乘（1）矩阵乘矩阵乘与数组乘有着较大差别，运算结果也完全不同。矩阵乘的运算符为“*”，运算是按矩阵的乘法规则进行，即参与乘运算的两矩阵的内维必须相同。设A、B为参与乘运算的 =A m×k B k×n。因此，参与运两矩阵，C为A和B的矩阵乘的结果，则它们必须满足关系C m ×n 算的两矩阵的顺序不能任意调换，因为A*B和B*A计算结果很可能是完全不一样的。如：A=[1 1 1;2 2 2;3 3 3]; B=A;

A*B ans= 6 6 6 12 12 12 18 18 18 F=ones(1,3); G=ones(3,1); F*G ans 3 G*F ans= 1 1 1 1 1 1 1 1 1 （2）数组乘数组乘的运算符为“.*”，运算符中的点号不能遗漏，也不能随意加空格符。参加数组乘运算的两数组的大小必须相等（即同维数组）。数组乘的结果是将两同维数组（矩阵）的对应元素逐一相乘，因此，A.*B和B.*A的计算结果是完全相同的，如： A=[1 1 1 1 1;2 2 2 2 2;3 3 3 3 3]; B=A; A.*B ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 B.*A ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 由于矩阵运算和数组运算的差异，能进行数组乘运算的两矩阵，不一定能进行矩阵乘运算。如 A=ones(1,3); B=A; A.*B ans= 1 1 1 A*A ???Error using= =>

数值分析matlab代码

1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6) disp('牛顿法'); i=1; n0=180; p0=pi/3; tol=10^(-6); for i=1:n0 p=p0-(p0-sin(p0))/(1-cos(p0)); if abs(p-p0)<=10^(-6) disp('用牛顿法求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次牛顿迭代后无法求出方程的解') end 2、disp('Steffensen加速'); p0=pi/3; for i=1:n0 p1=0.5*p0+0.5*cos(p0); p2=0.5*p1+0.5*cos(p1); p=p0-((p1-p0).^2)./(p2-2.*p1+p0); if abs(p-p0)<=10^(-6) disp('用Steffensen加速求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次Steffensen加速后无法求出方程的解') end 1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根， %误差限为 e=10^-4 disp('二分法')

a=0.2;b=0.26; tol=0.0001; n0=10; fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1; for i=1:n0 p=(a+b)/2; fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1; if fp==0||(abs((b-a)/2)0 a=p; else b=p; end end if i==n0&&~(fp==0||(abs((b-a)/2)

数值分析幂法与反幂法-matlab程序

数值分析幂法与反幂法 matlab程序随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量。要求 1）比较不同的原点位移和初值说明收敛性 2）给出迭代结果，生成DOC文件。 3）程序清单，生成M文件。解答： >> A=rand(5) %随机产生5*5矩阵求随机矩阵 A = 0.7094 0.1626 0.5853 0.6991 0.1493 0.7547 0.1190 0.2238 0.8909 0.2575 0.2760 0.4984 0.7513 0.9593 0.8407 0.6797 0.9597 0.2551 0.5472 0.2543 0.6551 0.3404 0.5060 0.1386 0.8143 >> B=A+A' %A矩阵和A的转置相加，得到随机对称矩阵B B = 1.4187 0.9173 0.8613 1.3788 0.8044 0.9173 0.2380 0.7222 1.8506 0.5979 0.8613 0.7222 1.5025 1.2144 1.3467 1.3788 1.8506 1.2144 1.0944 0.3929 0.8044 0.5979 1.3467 0.3929 1.6286

B=?? ????? ???? ?? ???6286.13929.03467.15979.08044 .03929.00944 .12144.18506 .13788.13467.12144.15025.17222.08613.05979.08506.17222.02380.09173.08044.03788.18613 .09173 .04187.1 编写幂法、反幂法程序： function [m,u,index,k]=pow(A,u,ep,it_max) % 求矩阵最大特征值的幂法，其中 % A 为矩阵； % ep 为精度要求，缺省为1e-5； % it_max 为最大迭代次数，缺省为100； % m 为绝对值最大的特征值； % u 为对应最大特征值的特征向量； % index ，当index=1时，迭代成功，当index=0时，迭代失败 if nargin<4 it_max=100; end if nargin<3 ep=1e-5; end n=length(A); index=0; k=0; m1=0; m0=0.01; % 修改移位参数，原点移位法加速收敛，为0时，即为幂法 I=eye(n) T=A-m0*I while k<=it_max v=T*u; [vmax,i]=max(abs(v)); m=v(i); u=v/m; if abs(m-m1)

MATLAB与数值分析实验报告一

MATLAB与数值分析实验报告报告人：秦旸照学号： 2015020901033 时间： 2016.4.8 电子科技大学电子工程学院

一、实验目的实验一：MATLAB软件平台与程序设计实验二、实验原理 1.熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算操作。（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序） 2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。（用.m文件编写进行符号因式分解和函数求反的程序） 3. 掌握Matlab函数的编写规范。 4.掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、三维曲线和面的填充、三维等高线等。（用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释） 5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。三、实验方案 1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。并以 x，y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图，允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。不允许使用sort函数。四、实验结果 1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。并以 x，y为坐标显示图像

数值分析算法在matlab中的实现

数值分析matlab实现高斯消元法： function[RA,RB,n,X]=gaus(A,b) B=[A b];n=length(b);RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1 for k=p+1:n m=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end else disp('请注意：因为RA=RB0, disp('请注意：因为RA~=RB，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1

matlab数值分析例题

1、在MATLAB 中用Jacobi 迭代法讨论线性方程组， 1231231234748212515 x x x x x x x x x -+=?? -+=-??-++=? （1）给出Jacobi 迭代法的迭代方程，并判定Jacobi 迭代法求解此方程组是否收敛。（2）若收敛，编程求解该线性方程组。解（1）：A=[4 -1 1;4 -8 1;-2 1 5] %线性方程组系数矩阵 A = 4 -1 1 4 -8 1 -2 1 5 >> D=diag(diag(A)) D = 4 0 0 0 -8 0 0 0 5 >> L=-tril(A,-1) % A 的下三角矩阵 L = 0 0 0 -4 0 0 2 -1 0 >> U=-triu(A,1) % A 的上三角矩阵 U = 0 1 -1 0 0 -1 0 0 0 B=inv(D)*(L+U) % B 为雅可比迭代矩阵 B = 0 0.2500 -0.2500 0.5000 0 0.1250 0.4000 -0.2000 0 >> r=eigs(B,1) %B 的谱半径

r = 0.3347 < 1 Jacobi迭代法收敛。 (2)在matlab上编写程序如下： A=[4 -1 1;4 -8 1;-2 1 5]; >> b=[7 -21 15]'; >> x0=[0 0 0]'; >> [x,k]=jacobi(A,b,x0,1e-7) x = 2.0000 4.0000 3.0000 k = 17 附jacobi迭代法的matlab程序如下： function [x,k]=jacobi(A,b,x0,eps) % 采用Jacobi迭代法求Ax=b的解 % A为系数矩阵 % b为常数向量 % x0为迭代初始向量 % eps为解的精度控制 max1= 300; %默认最多迭代300，超过300次给出警告D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\(L+U); f=D\b; x=B*x0+f; k=1; %迭代次数 while norm(x-x0)>=eps x0=x; x=B*x0+f; k=k+1; if(k>=max1) disp('迭代超过300次，方程组可能不收敛'); return; end end

数值分析matlab函数资料

1．求数值积分： fx=@(x)exp(1./x); >> quadl(fx,1,5) 2．获取x=xlsread('oillack.xls','sheet1','a1:a73') excel文件名是oillack.xls，sheet1是表名，a1:a73'是a列的1到73行 long x=xlsread('F:\学习\大三\大三下\巷道力学模型\新建文件夹(2)\1.xlsx','sheet1','a2:a') 3. 在matlab的图中插入文本框后将文本框旋转的方法： text(0.5,0.6,'渗透率/mD','Rotation',90) 4. matlab中插入一条直线的方法： line([0.01 0.01],[0 1.75]) 5．Matlab 中画三维图 x=-7.5:0.5:7.5; y=x; % 先产生x及y二个阵列 >> [x,y]=meshgrid(x,y); % 再以meshgrid形成二维的网格数据 >> z=x.^2+y.^2; % 产生z轴的数据 >> mesh(x,y,z) % 将z轴的变化值以网格方式画出 >> surf(x,y,z) % 将z轴的变化值以曲面方式画出 Matlab指数拟合方法 x=[1982 1992 2002]； y=[103.5 34.5 23.3]； cftool（x，y）在弹出的对话框选择fitting，弹出新的对话框选择new fit，然后在第三个下拉菜单（Type of fit）中选择Exponential，然后点击Apply，即可；最后结果 General model Exp1: f(x) = a*exp(b*x) Coefficients (with 95% confidence bounds): a = 1.453e+082 (-7.288e+084, 7.317e+084) b = -0.09312 (-0.3464, 0.1602)

数值分析实验— MATLAB实现

数值分析实验 ——MATLAB实现姓名sumnat 学号2013326600000 班级13级应用数学2班指导老师 2016年1月

一、插值：拉格朗日插值 (1) 1、代码： (1) 2、示例： (1) 二、函数逼近：最佳平方逼近 (2) 1、代码： (2) 2、示例： (2) 三、数值积分：非反常积分的Romberg算法 (3) 1、代码： (3) 2、示例： (4) 四、数值微分：5点法 (5) 1、代码： (5) 2、示例： (6) 五、常微分方程：四阶龙格库塔及Adams加速法 (6) 1、代码：四阶龙格库塔 (6) 2、示例： (7) 3、代码：Adams加速法 (7) 4、示例： (8) 六、方程求根：Aitken 迭代 (8) 1、代码： (8) 2、示例： (9) 七、线性方程组直接法：三角分解 (9) 1、代码： (9) 2、示例： (10) 八、线性方程组迭代法：Jacobi法及G-S法 (11) 1、代码：Jacobi法 (11) 2、示例： (12) 3、代码：G-S法 (12) 4、示例： (12) 九、矩阵的特征值及特征向量：幂法 (13) 1、代码： (13) 2、示例： (13)

一、插值：拉格朗日插值 1、代码： function z=LGIP(x,y)%拉格朗日插值 n=size(x); n=n(2);%计算点的个数 syms a; u=0;%拉格朗日多项式 f=1;%插值基函数 for i=1:n for j=1:n if j==i f=f; else f=f*(a-x(j))/(x(i)-x(j)); end end u=u+y(i)*f;f=1; end z=expand(u);%展开 2、示例： >> x=1:6; y1=x.^5+3*x.^2-6; y2=sin(x)+sqrt(x); >> f1=LGIP(x,y1) f1 = -6+3*a^2+a^5 %可知多项式吻合得很好 >> f2=vpa(LGIP(x,y2),3) f2 = .962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5

数理方程基于matlab的数值解法

数理方程数值解法与其在matlab软件上的实现张体强1026222 廖荣发1026226 [摘要] 数学物理方程的数值解在实际生活中越来越使用，首先基于偏微分数值解的思想上，通过matlab软件的功能，研究其数学物理方程的数值解，并通过对精确解和数值解进行对比，追究其数值解的可行性，在此，给出相关例子和程序代码，利于以后的再次研究和直接使用。 [关键字] 偏微分方程数值解matlab 数学物理方程的可视化一：研究意义在我们解数学物理方程，理论上求数学物理方程的定解有着多种解法，但是有许多定解问题却不能严格求解，只能用数值方法求出满足实际需要的近似解。而且实际问题往往很复杂，这时即便要解出精确解就很困难，有时甚至不可能，另一方面，在建立数学模型时，我们已作了很多近似，所以求出的精确解也知识推导出的数学问题的精确解，并非真正实际问题的精确解。因此，我们有必要研究近似解法，只要使所求得的近似解与精确解之间的误差在规定的范围内，则仍能满足实际的需要，有限差分法和有限元法是两种最常用的

求解数学物理方程的数值解法,而MATLAB 在这一方面具有超强的数学功能，可以用来求其解。二：数值解法思想和步骤 2.1:网格剖分为了用差分方法求解上述问题，将求解区域 {}(,)|01,01x t x t Ω=≤≤≤≤作剖分。将空间区间[0,1]作m 等分，将时间[0,1]区间作n 等分，并记 1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。分别称h 和τ 为空间和时间步长。用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。 2.2：差分格式的建立 0u u t x ??-=??………………………………（1）设G 是,x t 平面任一有界域，据Green 公式(参考数学物理方程第五章)： ( )()G u u dxdt udt udx t x Γ??-=--??? ? 其中G Γ=?。于是可将（1）式写成积分守恒形式： ()0udt udx Γ --=? (2) 我们先从（2）式出发构造熟知的Lax 格式设网格如下图所示

数值分析五个题目的C语言及Matlab程序

本文档包含上一个文档中的五个数值分析实验题C语言程序及Matlab程序实验一 C程序 #include "stdio.h" #include "math.h" void main() { inti=0; float a=0.1,b=1.9,t=0.0,e=1.9; if((pow(a,7)-28*pow(a,4)+14)*(pow(b,7)-28*pow(b,4)+14)<0) if((7*pow(x,6)-112*pow(x,3))) printf("x=%f,i=%d,e=%f\n",x,i,e); for(i=1;i<7&&e>0.00001;i++) { t=x; x=x-(pow(x,7)-28*pow(x,4)+14)/(7*pow(x,6)-112*pow(x,3)); e=fabs(t-x); printf("x=%f,i=%d,e=%f\n",x,i,e);

} } Matable 程序 i=0; x=1.9;t=0.0;e=1.9; disp(['i=',num2str(i),' ','x=',num2str(x),' ','e=',num2str(e)]); for i=1:7 t=x; x=x-(x^7-28*x^4+14)/(7*x^6-112*x^3); e=abs(t-x); disp(['i=',num2str(i),' ','x=',num2str(x),' ','e=',num2str(e)]); if e<0.00001 break; end end 实验二 C程序 #include"stdio.h" #include"math.h" //已知量 double x[10]={1,2,3,4,5,6,7,8,9,10};