3．钢水包使用次数多以后，钢包的容积增大，数据如下：
| x | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
| y | 6.42 | 8.2 | 9.58 | 9.5 | 9.7 | 10 | 9.93 | 9.99 |
| 10 | 11 | 12 | 13 | 14 | 15 | 16 |
| 10.49 | 10.59 | 10.60 | 10.8 | 10.6 | 10.9 | 10.76 |
试从中找出使用次数和容积之间的关系，计算均方差。（用
拟合）
- 为什么要改提指标？
- 如何去做？
因为用经验公式拟合以后得到的方程应该是一个线性的方程，之后求解其驻点。
但是现在的公式是非线性的。
所以需要把分母移到y的一侧，再移到等式右边，将等式右边x与y的系数求和， 记为新的方程组，对这个方程组求解他们的参数的偏导数。
方法一
x0=2:16;
y0=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.8 10.6 10.9 10.76];
x0=x0';
y0=y0';
yx=@(p,x)p(1)+(p(2)-p(1).*p(3))./(x+p(3));
p=lsqcurvefit(yx,rand(3,1),x0,y0);
fprintf('equation is %0.2f+(%0.2f-%0.2f*%0.2f)/(x+%0.2f)',p(1),p(2),p(1),p(3),p(3))
y=p(1)+(p(2)-p(1).*p(3))./(x0+p(3));
plot(x0,y);
hold on
plot(x0,y0,'xg')
var=0;
for i=1:15.
var=var+(y0(i)-(p(1)+(p(2)-p(1).*p(3))./(x0(i)+p(3))))^2;
end
var=sqrt(var/15)
equation is 11.22+(-12.22-11.22*-0.41)/(x+-0.41)
var =
0.1938
方法二：
将分母移到y一侧得到：
令
对f求a, b, c的偏导后等于0，化简有：
根据矩阵:
A为系数矩阵。
因此，inv(A)*B可以求解a, b, c
x=2:16;
y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.8 10.6 10.9 10.76];
A=[sum(x.^2) sum(x) -sum(x.*y);sum(x) 15 -sum(y); sum(x.*y) sum(y) -sum(y.^2)];
B=[sum(x.^2.*y) ;sum(x.*y); sum(y.^2.*x)];
theta=inv(A)*B;
n=2:0.01:18;
a=theta(1);
b=theta(2);
c=theta(3);
y_hat=(a*n+b)./(n+c);
plot(n,y_hat)
hold on
plot(x,y,'xg');
var=0;
for i=1:15
var=var+(y(i)-(a*x(i)+b)./(x(i)+c))^2;
end
var=sqrt(var/15);
fprintf('equation is y = (%0.2fx+%0.2f)/(x+%0.2f)\n',a,b,c)
fprintf('the var is %0.2f\n',var)
equation is y = (11.34x+-12.50)/(x+-0.34)
the var is 0.22
1、
A=[4 -1 0 -1 0 0;
-1 4 -1 0 -1 0;
0 -1 4 -1 0 -1;
-1 0 -1 4 -1 0;
0 -1 0 -1 4 -1;
0 0 -1 0 -1 4;];
b=[0 5 -2 5 -2 6]';
x0=[0 0 0 0 0 0]';
err=10^(-4);
D = diag(diag(A));
B = D\(D-A);
p=vrho(B)
x1=[1 1 1 1 1 1]';
c=(x0-x1).^2;
d=sqrt(sum(c(:)));
i=0;
while d>err
i=i+1;
x1(1)=(x0(2)+x0(4))/4;
x1(2)=(5+x0(1)+x0(3)+x0(5))/4;
x1(3)=(x0(2)+x0(4)+x0(6)-2)/4;
x1(4)=(5+x0(1)+x0(3)+x0(5))/4;
x1(5)=(x0(2)+x0(4)+x0(6)-2)/4;
x1(6)=(6+x0(3)+x0(5))/4;
c=(x0-x1).^2;
d=sqrt(sum(c(:)));
x0=x1;
end
fprintf('Jacobi Function:\n')
fprintf('approximation is x1=%0.2f\nx2=%0.2f\nx3=%0.2f\nx4=%0.2f\nx5=%0.2f\nx6=%0.2f\n',x1(1),x1(2),x1(3),x1(4),x1(5),x1(6));
fprintf('iteration is %d\n',i);
p =
0.6830
Jacobi Function:
approximation is x1=1.00
x2=2.00
x3=1.00
x4=2.00
x5=1.00
x6=2.00
iteration is 28
谱半径为0.683<1 验证了后续的收敛。
2、
B=inv(tril(A,0))*triu(A,1);
p=vrho(B)
x1=[1 1 1 1 1 1]';
c=(x0-x1).^2;
d=sqrt(sum(c(:)));
i=0;
while d>err
i=i+1;
x1(1)=(x0(2)+x0(4))/4;
x1(2)=(5+x1(1)+x0(3)+x0(5))/4;
x1(3)=(x1(2)+x0(4)+x0(6)-2)/4;
x1(4)=(5+x1(1)+x1(3)+x0(5))/4;
x1(5)=(x1(2)+x1(4)+x0(6)-2)/4;
x1(6)=(6+x1(3)+x1(5))/4;
c=(x0-x1).^2;
d=sqrt(sum(c(:)));
x0=x1;
end
fprintf('Gauss-Seidel Function:\n')
fprintf('approximation is x1=%0.2f\nx2=%0.2f\nx3=%0.2f\nx4=%0.2f\nx5=%0.2f\nx6=%0.2f\n',x1(1),x1(2),x1(3),x1(4),x1(5),x1(6));
fprintf('iteration is %d\n',i);
p =
0.4806
Gauss-Seidel Function:
approximation is x1=1.00
x2=2.00
x3=1.00
x4=2.00
x5=1.00
x6=2.00
iteration is 15
谱半径为0.4806<1 验证了后续的收敛。
3、
A=[4 -1 0 -1 0 0;
-1 4 -1 0 -1 0;
0 -1 4 -1 0 -1;
-1 0 -1 4 -1 0;
0 -1 0 -1 4 -1;
0 0 -1 0 -1 4;];
b=[0 5 -2 5 -2 6]';
x0=[0 0 0 0 0 0]';
err=10^(-4);
x1=[1 1 1 1 1 1]';
c=(x0-x1).^2;
d=sqrt(sum(c(:)));
w=0;
h=0.1;
step=1;
fprintf('SOR Function:\n')
while 1
w=w+h;
if w>=2;
break
end
i=0;
while d>err
i=i+1;
x1(1)=x0(1)+w*(x0(2)+x0(4)-4*x0(1))/4;
x1(2)=x0(2)+w*(5+x1(1)+x0(3)+x0(5)-4*x0(2))/4;
x1(3)=x0(3)+w*(x1(2)+x0(4)+x0(6)-2-4*x0(3))/4;
x1(4)=x0(4)+w*(5+x1(1)+x1(3)+x0(5)-4*x0(4))/4;
x1(5)=x0(5)+w*(x1(2)+x1(4)+x0(6)-2-4*x0(5))/4;
x1(6)=x0(6)+w*(6+x1(3)+x1(5)-4*x0(6))/4;
c=(x0-x1).^2;
d=sqrt(sum(c(:)));
x0=x1;
end
w_sum(step)=i;
step=step+1;
fprintf('iteration is %d ;w is %0.1f\n',i,w);
fprintf('approximation is x1=%0.2f x2=%0.2f x3=%0.2f x4=%0.2f x5=%0.2f x6=%0.2f\n',x1(1),x1(2),x1(3),x1(4),x1(5),x1(6));
x0=[0 0 0 0 0 0]';
x1=[1 1 1 1 1 1]';
c=(x0-x1).^2;
d=sqrt(sum(c(:)));
end
[min_value,min_ind]=min(w_sum);
fprintf('the minimum w is %0.1f, its value is %d\n',min_ind/10,min_value)
SOR Function:
iteration is 214 ;w is 0.1
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 112 ;w is 0.2
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 75 ;w is 0.3
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 55 ;w is 0.4
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 43 ;w is 0.5
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 34 ;w is 0.6
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 28 ;w is 0.7
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 23 ;w is 0.8
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 19 ;w is 0.9
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 15 ;w is 1.0
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 12 ;w is 1.1
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 10 ;w is 1.2
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 12 ;w is 1.3
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 16 ;w is 1.4
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 20 ;w is 1.5
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 26 ;w is 1.6
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 36 ;w is 1.7
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 57 ;w is 1.8
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
iteration is 117 ;w is 1.9
approximation is x1=1.00 x2=2.00 x3=1.00 x4=2.00 x5=1.00 x6=2.00
the minimum w is 1.2, its value is 10
方程组一：
function [k1,k2]= T6K(y1,y2,x)
k1=-2*y1+y2+2*sin(x);
k2=y1-2*y2+2*cos(x)-2*sin(x);
clear;clc;
y1(1)=2;
y2(1)=3;
h=0.1;
x(1)=0;
n=0:0.1:10;
Y1=2*exp(-n)+sin(n);
Y2=2*exp(-n)+cos(n);
for i=1:100
[k11,k12]=T6K(y1(i),y2(i),x(i));
[k21,k22]=T6K((y1(i)+k11*h),(y2(i)+k12*h),(x(i)+h));
y1(i+1)=y1(i)+h/2*(k11+k21);
y2(i+1)=y2(i)+h/2*(k12+k22);
x(i+1)=x(i)+h;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
end
i=i+1;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
RW=[x(1),Y1(1),y1(1),err1(1)];
for i=1:10/h
if mod(i,10)==0
RW=[RW;x(i),Y1(i),y1(i),err1(i)];
end
end
fprintf('\t x\texact y1 approximate y1 err y1\n')
disp(RW)
RW=[x(1),Y2(1),y2(1),err2(1)];
for i=1:10/h
if mod(i,10)==0
RW=[RW;x(i),Y2(i),y2(i),err2(i)];
end
end
fprintf('\n\n\t x\texact y2 approximate y2 err y2\n')
disp(RW)
x exact y1 approximate y1 err y1
0 2.0000 2.0000 0
0.9000 1.5965 1.5968 0.0003
1.9000 1.2454 1.2439 0.0016
2.9000 0.3493 0.3481 0.0012
3.9000 -0.6473 -0.6463 0.0010
4.9000 -0.9676 -0.9649 0.0027
5.9000 -0.3684 -0.3663 0.0021
6.9000 0.5805 0.5802 0.0003
7.9000 0.9997 0.9973 0.0024
8.9000 0.5013 0.4990 0.0023
9.9000 -0.4574 -0.4575 0.0001
x exact y2 approximate y2 err y2
0 3.0000 3.0000 0
0.9000 1.4347 1.4348 0.0001
1.9000 -0.0242 -0.0226 0.0015
2.9000 -0.8609 -0.8583 0.0026
3.9000 -0.6854 -0.6834 0.0020
4.9000 0.2014 0.2014 0.0000
5.9000 0.9330 0.9312 0.0018
6.9000 0.8177 0.8159 0.0019
7.9000 -0.0453 -0.0454 0.0002
8.9000 -0.8652 -0.8635 0.0017
9.9000 -0.8891 -0.8871 0.0020
步长h=0.001
clear;clc;
y1(1)=2;
y2(1)=3;
h=0.001;
x(1)=0;
n=0:h:10;
Y1=2*exp(-n)+sin(n);
Y2=2*exp(-n)+cos(n);
for i=1:10/h
[k11,k12]=T6K(y1(i),y2(i),x(i));
[k21,k22]=T6K((y1(i)+k11*h),(y2(i)+k12*h),(x(i)+h));
y1(i+1)=y1(i)+h/2*(k11+k21);
y2(i+1)=y2(i)+h/2*(k12+k22);
x(i+1)=x(i)+h;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
end
i=i+1;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
RW=[x(1),Y1(1),y1(1),err1(1)];
for i=1:10/h
if mod(i,500)==0
RW=[RW;x(i),Y1(i),y1(i),err1(i)];
end
end
fprintf('\t x\texact y1 approximate y1 err y1\n')
disp(RW)
RW=[x(1),Y2(1),y2(1),err2(1)];
for i=1:10/h
if mod(i,500)==0
RW=[RW;x(i),Y2(i),y2(i),err2(i)];
end
end
fprintf('\n\n\t x\texact y2 approximate y2 err y2\n')
disp(RW)
x exact y1 approximate y1 err y1
0 2.0000 2.0000 0
0.4990 1.6928 1.6928 0.0000
0.9990 1.5774 1.5774 0.0000
1.4990 1.4441 1.4441 0.0000
1.9990 1.1807 1.1807 0.0000
2.4990 0.7636 0.7636 0.0000
2.9990 0.2418 0.2418 0.0000
3.4990 -0.2894 -0.2894 0.0000
3.9990 -0.7195 -0.7195 0.0000
4.4990 -0.9551 -0.9551 0.0000
4.9990 -0.9457 -0.9457 0.0000
5.4990 -0.6981 -0.6981 0.0000
5.9990 -0.2754 -0.2754 0.0000
6.4990 0.2172 0.2172 0.0000
6.9990 0.6581 0.6581 0.0000
7.4990 0.9388 0.9388 0.0000
7.9990 0.9902 0.9902 0.0000
8.4990 0.7995 0.7995 0.0000
8.9990 0.4133 0.4133 0.0000
9.4990 -0.0740 -0.0740 0.0000
9.9990 -0.5431 -0.5431 0.0000
x exact y2 approximate y2 err y2
0 3.0000 3.0000 0
0.4990 2.0923 2.0923 0.0000
0.9990 1.2776 1.2776 0.0000
1.4990 0.5184 0.5184 0.0000
1.9990 -0.1443 -0.1443 0.0000
2.4990 -0.6362 -0.6362 0.0000
2.9990 -0.8902 -0.8902 0.0000
3.4990 -0.8764 -0.8764 0.0000
3.9990 -0.6177 -0.6177 0.0000
4.4990 -0.1895 -0.1895 0.0000
4.9990 0.2962 0.2962 0.0000
5.4990 0.7161 0.7161 0.0000
5.9990 0.9649 0.9649 0.0000
6.4990 0.9798 0.9798 0.0000
6.9990 0.7564 0.7564 0.0000
7.4990 0.3487 0.3487 0.0000
7.9990 -0.1438 -0.1438 0.0000
8.4990 -0.6008 -0.6008 0.0000
8.9990 -0.9105 -0.9105 0.0000
9.4990 -0.9971 -0.9971 0.0000
9.9990 -0.8395 -0.8395 0.0000
（2）
function [k1,k2]= T6Ky2(y1,y2,x)
k1=-2*y1+y2+2*sin(x);
k2=998*y1-999*y2+999*cos(x)-999*sin(x);
欧拉显式的方法：
clear;clc;
y1(1)=2;
y2(1)=3;
h=0.1;
x(1)=0;
n=0:h:10;
Y1=2*exp(-n)+sin(n);
Y2=2*exp(-n)+cos(n);
for i=1:10/h
[k11,k12]=T6Ky2(y1(i),y2(i),x(i));
[k21,k22]=T6Ky2((y1(i)+k11*h),(y2(i)+k12*h),(x(i)+h));
y1(i+1)=y1(i)+h/2*(k11+k21);
y2(i+1)=y2(i)+h/2*(k12+k22);
x(i+1)=x(i)+h;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
end
i=i+1;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
x=single(x);
RW=[x(1),Y1(1),y1(1),err1(1)];
for i=1:10/h
if mod(i,10)==0
RW=[RW;x(i),Y1(i),y1(i),err1(i)];
end
end
fprintf('\t x\texact y1 approximate y1 err y1\n')
disp(RW)
RW=[x(1),Y2(1),y2(1),err2(1)];
for i=1:10/h
if mod(i,10)==0
RW=[RW;x(i),Y2(i),y2(i),err2(i)];
end
end
fprintf('\n\n\t x\texact y2 approximate y2 err y2\n')
disp(RW)
x exact y1 approximate y1 err y1
1.0e+25 *
0 0.0000 0.0000 0
0.0000 0.0000 8.0442 8.0442
0.0000 0.0000 Inf Inf
0.0000 0.0000 Inf Inf
0.0000 -0.0000 Inf Inf
0.0000 -0.0000 Inf Inf
0.0000 -0.0000 Inf Inf
0.0000 0.0000 Inf Inf
0.0000 0.0000 Inf Inf
0.0000 0.0000 NaN NaN
0.0000 -0.0000 NaN NaN
x exact y2 approximate y2 err y2
1.0e+28 *
0 0.0000 0.0000 0
0.0000 0.0000 -8.0281 8.0281
0.0000 -0.0000 -Inf Inf
0.0000 -0.0000 -Inf Inf
0.0000 -0.0000 -Inf Inf
0.0000 0.0000 -Inf Inf
0.0000 0.0000 -Inf Inf
0.0000 0.0000 -Inf Inf
0.0000 -0.0000 -Inf Inf
0.0000 -0.0000 NaN NaN
0.0000 -0.0000 NaN NaN
clear;clc;
y1(1)=2;
y2(1)=3;
h=0.001;
x(1)=0;
n=0:h:10;
Y1=2*exp(-n)+sin(n);
Y2=2*exp(-n)+cos(n);
for i=1:10/h
[k11,k12]=T6Ky2(y1(i),y2(i),x(i));
[k21,k22]=T6Ky2((y1(i)+k11*h),(y2(i)+k12*h),(x(i)+h));
y1(i+1)=y1(i)+h/2*(k11+k21);
y2(i+1)=y2(i)+h/2*(k12+k22);
x(i+1)=x(i)+h;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
end
i=i+1;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
x=single(x);
RW=[x(1),Y1(1),y1(1),err1(1)];
for i=1:10/h
if mod(i,500)==0
RW=[RW;x(i),Y1(i),y1(i),err1(i)];
end
end
fprintf('\t x\texact y1 approximate y1 err y1\n')
disp(RW)
RW=[x(1),Y2(1),y2(1),err2(1)];
for i=1:10/h
if mod(i,500)==0
RW=[RW;x(i),Y2(i),y2(i),err2(i)];
end
end
fprintf('\n\n\t x\texact y2 approximate y2 err y2\n')
disp(RW)
x exact y1 approximate y1 err y1
0 2.0000 2.0000 0
0.4990 1.6928 1.6928 0.0000
0.9990 1.5774 1.5774 0.0000
1.4990 1.4441 1.4441 0.0000
1.9990 1.1807 1.1807 0.0000
2.4990 0.7636 0.7636 0.0000
2.9990 0.2418 0.2418 0.0000
3.4990 -0.2894 -0.2894 0.0000
3.9990 -0.7195 -0.7195 0.0000
4.4990 -0.9551 -0.9551 0.0000
4.9990 -0.9457 -0.9457 0.0000
5.4990 -0.6981 -0.6981 0.0000
5.9990 -0.2754 -0.2754 0.0000
6.4990 0.2172 0.2172 0.0000
6.9990 0.6581 0.6581 0.0000
7.4990 0.9388 0.9388 0.0000
7.9990 0.9902 0.9902 0.0000
8.4990 0.7995 0.7995 0.0000
8.9990 0.4133 0.4133 0.0000
9.4990 -0.0740 -0.0740 0.0000
9.9990 -0.5431 -0.5431 0.0000
x exact y2 approximate y2 err y2
0 3.0000 3.0000 0
0.4990 2.0923 2.0923 0.0000
0.9990 1.2776 1.2776 0.0000
1.4990 0.5184 0.5184 0.0000
1.9990 -0.1443 -0.1443 0.0000
2.4990 -0.6362 -0.6362 0.0000
2.9990 -0.8902 -0.8902 0.0000
3.4990 -0.8764 -0.8764 0.0000
3.9990 -0.6177 -0.6177 0.0000
4.4990 -0.1895 -0.1895 0.0000
4.9990 0.2962 0.2962 0.0000
5.4990 0.7161 0.7161 0.0000
5.9990 0.9649 0.9649 0.0000
6.4990 0.9798 0.9798 0.0000
6.9990 0.7564 0.7564 0.0000
7.4990 0.3487 0.3487 0.0000
7.9990 -0.1438 -0.1438 0.0000
8.4990 -0.6008 -0.6008 0.0000
8.9990 -0.9105 -0.9105 0.0000
9.4990 -0.9971 -0.9971 0.0000
9.9990 -0.8395 -0.8395 0.0000
中点隐式：
clear;clc;
h = 0.1;
x(1) = 0;
y1(1)= 2;
y2(1)= 3;
x=0:h:10;
Y1=2*exp(-x)+sin(x);
Y2=2*exp(-x)+cos(x);
I=eye(2);
A=[-2 1;998 -999];
A1=I-h/2*A;
for i= 1: 10/h
[K1(1),K1(2)]=T6Ky2(y1(i),y2(i),x(i));
y_m1=y1(i)+h/2*K1(1);
y_m2=y2(i)+h/2*K1(2);
[D(1),D(2)]=T6Ky2(y_m1,y_m2,x(i));
K2=A1 \ D';
y1(i+1)=y1(i)+(K1(1)+K2(1))*h/2;
y2(i+1)=y2(i)+(K1(2)+K2(2))*h/2;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
end
i=i+1;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
RW=[x(1),Y1(1),y1(1),err1(1)];
for i=1:10/h
if mod(i,10)==0
RW=[RW;x(i),Y1(i),y1(i),err1(i)];
end
end
fprintf('\t x\texact y1 approximate y1 err y1\n')
disp(RW)
RW=[x(1),Y2(1),y2(1),err2(1)];
for i=1:10/h
if mod(i,10)==0
RW=[RW;x(i),Y2(i),y2(i),err2(i)];
end
end
fprintf('\n\n\t x\texact y2 approximate y2 err y2\n')
disp(RW)
x exact y1 approximate y1 err y1
0 2.0000 2.0000 0
0.9000 1.5965 1.5954 0.0010
1.9000 1.2454 1.2448 0.0007
2.9000 0.3493 0.3494 0.0001
3.9000 -0.6473 -0.6468 0.0005
4.9000 -0.9676 -0.9673 0.0003
5.9000 -0.3684 -0.3687 0.0003
6.9000 0.5805 0.5799 0.0006
7.9000 0.9997 0.9993 0.0004
8.9000 0.5013 0.5014 0.0002
9.9000 -0.4574 -0.4569 0.0006
x exact y2 approximate y2 err y2
0 3.0000 3.0000 0
0.9000 1.4347 1.4337 0.0010
1.9000 -0.0242 -0.0248 0.0007
2.9000 -0.8609 -0.8609 0.0001
3.9000 -0.6854 -0.6850 0.0005
4.9000 0.2014 0.2017 0.0003
5.9000 0.9330 0.9327 0.0003
6.9000 0.8177 0.8172 0.0006
7.9000 -0.0453 -0.0457 0.0004
8.9000 -0.8652 -0.8650 0.0001
9.9000 -0.8891 -0.8885 0.0006
clear;clc;
h = 0.001;
x(1) = 0;
y1(1)= 2;
y2(1)= 3;
x=0:h:10;
Y1=2*exp(-x)+sin(x);
Y2=2*exp(-x)+cos(x);
I=eye(2);
A=[-2 1;998 -999];
A1=I-h/2*A;
for i= 1: 10/h
[K1(1),K1(2)]=T6Ky2(y1(i),y2(i),x(i));
y_m1=y1(i)+h/2*K1(1);
y_m2=y2(i)+h/2*K1(2);
[D(1),D(2)]=T6Ky2(y_m1,y_m2,x(i));
K2=A1 \ D';
y1(i+1)=y1(i)+(K1(1)+K2(1))*h/2;
y2(i+1)=y2(i)+(K1(2)+K2(2))*h/2;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
end
i=i+1;
yreal1=2*exp(-x(i))+sin(x(i));
yreal2=2*exp(-x(i))+cos(x(i));
err1(i)=abs(y1(i)-yreal1);
err2(i)=abs(y2(i)-yreal2);
RW=[x(1),Y1(1),y1(1),err1(1)];
for i=1:10/h
if mod(i,500)==0
RW=[RW;x(i),Y1(i),y1(i),err1(i)];
end
end
fprintf('\t x\texact y1 approximate y1 err y1\n')
disp(RW)
RW=[x(1),Y2(1),y2(1),err2(1)];
for i=1:10/h
if mod(i,500)==0
RW=[RW;x(i),Y2(i),y2(i),err2(i)];
end
end
fprintf('\n\n\t x\texact y2 approximate y2 err y2\n')
disp(RW)
x exact y1 approximate y1 err y1
0 2.0000 2.0000 0
0.4990 1.6928 1.6928 0.0000
0.9990 1.5774 1.5774 0.0000
1.4990 1.4441 1.4441 0.0000
1.9990 1.1807 1.1807 0.0000
2.4990 0.7636 0.7636 0.0000
2.9990 0.2418 0.2418 0.0000
3.4990 -0.2894 -0.2894 0.0000
3.9990 -0.7195 -0.7195 0.0000
4.4990 -0.9551 -0.9551 0.0000
4.9990 -0.9457 -0.9457 0.0000
5.4990 -0.6981 -0.6981 0.0000
5.9990 -0.2754 -0.2754 0.0000
6.4990 0.2172 0.2172 0.0000
6.9990 0.6581 0.6581 0.0000
7.4990 0.9388 0.9388 0.0000
7.9990 0.9902 0.9902 0.0000
8.4990 0.7995 0.7995 0.0000
8.9990 0.4133 0.4133 0.0000
9.4990 -0.0740 -0.0740 0.0000
9.9990 -0.5431 -0.5431 0.0000
x exact y2 approximate y2 err y2
0 3.0000 3.0000 0
0.4990 2.0923 2.0923 0.0000
0.9990 1.2776 1.2776 0.0000
1.4990 0.5184 0.5184 0.0000
1.9990 -0.1443 -0.1443 0.0000
2.4990 -0.6362 -0.6362 0.0000
2.9990 -0.8902 -0.8902 0.0000
3.4990 -0.8764 -0.8764 0.0000
3.9990 -0.6177 -0.6177 0.0000
4.4990 -0.1895 -0.1895 0.0000
4.9990 0.2962 0.2962 0.0000
5.4990 0.7161 0.7161 0.0000
5.9990 0.9649 0.9649 0.0000
6.4990 0.9798 0.9798 0.0000
6.9990 0.7564 0.7564 0.0000
7.4990 0.3487 0.3487 0.0000
7.9990 -0.1438 -0.1438 0.0000
8.4990 -0.6008 -0.6008 0.0000
8.9990 -0.9105 -0.9105 0.0000
9.4990 -0.9971 -0.9971 0.0000
9.9990 -0.8395 -0.8395 0.0000
7.用有限差分法求解边值问题（h=0.1），并图形显示。
.
NOTE:
h1=0.1;
h2=0.0001;
y1=T7run(h1)
y2=T7run(h2);
n1=-1:h1:1;
n2=-1:h2:1;
plot(n1,y1,n2,y2)
legend('h=0.1','h=0.0001')
title('T7')
clear
function y=T7run(h)
% The whole coordinate axis is shifted to the right by two.
% y(-1)->y(1);y(1)->y(3)
y(1)=1;
y(3)=1;
len=2/h-1;
% create the diag
for i=1:len
x(i)=-1+i*h;
q(i)=(1+x(i)^2);
d(i)=q(i)+2/(h^2);
end
% create the d's sides
for i=1:len-1
d1(i)=-1/h^2;
end
% create the matrix A
A=diag(d1,1)+diag(d1,-1)+diag(d);
B(1)=y(1)/h^2;
B(len)=y(3)/h^2;
% use LU function when the h not enough small.
% [L,U]=lu(A);
% Ly=d
beta(1)=(-1/h^2)/A(1,1);
for i=2:len-1
beta(i)=(-1/h^2)/(A(i,i)-(-1/h^2)*beta(i-1));
end
y_m(1)=B(1)/A(1,1);
for i=2:len
y_m(i)=(B(i)+1/h^2*y_m(i-1))/(A(i,i)+(1/h^2)*beta(i-1));
end
% Ux=y
x1(len)=y_m(len);
for i=len-1:-1:1
x1(i)=y_m(i)-beta(i)*x1(i+1);
end
% y_m=inv(L)*B';
% x1=inv(U)*y_m;
y=[1,x1,1];
6091
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
