利用迭代法求解定非线性方程及方程组，使得误差不超过10^(-8)。同时应用迭代加速技术，提交迭代运算效率。
此题需要用到的MATLAB代码及附录：
附录6 二分法作根的隔离
% 附录6 二分法作根的隔离
%% 二分法作根的隔离
clear %清除变量
clc %清除命令行窗口代码
format long
aa=input('\n请输入自变量x的区间：\n');
x=[aa(1):0.1:aa(2)];
y=fun(x);
plot(x,y);
hold on
ezplot('0');
xlim([aa(1),aa(2)]);
k=0 ; %二分法迭代次数
e=aa(2)-aa(1); %区间长度大小
a=aa(1);b=aa(2);c=(a+b)/2;
err1=1e-4; err2=1e-4;
while e>=err1 && abs(fun(c))>=err2
if fun(a)*fun(c)<0
b=c;
end
if fun(a)*fun(c)>0
a=c;
end
if fun(a)*fun(c)==0
a=c;b=c;
end
c=(a+b)/2;
e=e/2;
k=k+1;
end
xx=c, %二分法求得的近似解
fid=fopen('data_xx1.txt','w');
fprintf(fid,'%.4f',xx); %保存二分法求得的近似解（满足1e-4精度要求）
fclose(fid);
fid=fopen('data_xx2.txt','w');
fprintf(fid,'%.4f',b); %保存二分法求得的近似解（不满足精度要求）用于两点弦割法
fclose(fid);
%% 应创建一个“fun.m”函数文件与第3节运行程序处于同一文件夹
function y = fun(x) %定义目标函数
y = x.^6-5*x.^5+3*x.^4+x.^3-7*x.^2+7*x-20;
end附录7 利用简单迭代法求解方程
% 附录7 利用简单迭代法求解方程
%% 利用简单迭代法求解方程
clear %清除变量
clc %清除命令行窗口代码
format long
err1=input('\n请输入最大允许误差：\n');
k=0 ; %迭代次数
x(1)=load('data_xx1.txt');
e=1.0+err1;
while abs(e)>=err1
k=k+1;
x(k+1)=(5*x(k)^5-3*x(k)^4-x(k)^3+7*x(k)^2-7*x(k)+20)^(1/6);
e=x(k+1)-x(k);
end
xx=x(k+1) %简单迭代法求得的近似解
fid=fopen('data_xx_results.txt','w');
fprintf(fid,'%.8f',xx);
fclose(fid);
附录8 利用牛顿法求解方程
% 附录8 利用牛顿法求解方程
%% 利用牛顿法求解方程
clear %清除变量
clc %清除命令行窗口代码
format long
err1=input('\n请输入最大允许误差：\n');
x(1)=load('data_xx1.txt');
e=1.0+err1;
syms t
fun1 =@(t) t.^6-5*t.^5+3*t.^4+t.^3-7*t.^2+7*t-20;
fun2 = diff(fun1(t));
fun2 = matlabFunction(fun2);
max1=1000; %最大迭代次数
y(1)=fun1(x(1));
for k = 1: max1
x(k+1) = x(k)-fun1(x(k))/fun2(x(k)); %Newton法的迭代公式
y(k+1) = fun1(x(k+1));
e = abs(x(k+1)-x(k)); %误差区间
if e<err1
break
end
end
xx=x(k+1) %牛顿迭代法求得的近似解
figure;
plot([1:k+1],y)
xlabel('X'); %横轴为迭代次数
ylabel('Y'); %纵轴为方程的值f(x)
grid;附录9 利用两点弦割法求解方程
% 附录9 利用两点弦割法求解方程
%% 利用两点弦割法求解方程
clear %清除变量
clc %清除命令行窗口代码
format long
err1 =input('\n请输入最大允许误差：\n');
x(1)=load('data_xx.txt'); y(1) = fun(x(1));
x(2)=x(1)+1e-5; y(2) = fun(x(2));
max1=1000; %最大迭代次数
for k = 2: max1
x(k+1) = x(k)-fun(x(k))*(x(k)-x(k-1))/(fun(x(k))-fun(x(k-1))); %两点弦割法的迭代公式
y(k+1) = fun(x(k+1));
e = abs(x(k+1)-x(k)); %误差区间
if e < err1
break
end
end
xx=x(k+1) %牛顿迭代法求得的近似解
figure;
plot([1:k+1],y)
xlabel('迭代步数');
ylabel('误差');
grid;
附录10 利用单点弦割法求解方程
% 附录10 利用单点弦割法求解方程
%% 利用单点弦割法求解方程
clear %清除变量
clc %清除命令行窗口代码
format long
err1 =input('\n请输入最大允许误差：\n');
x(1)=load('data_xx.txt'); y(1) = fun(x(1));
x(2)=x(1)+1e-2; y(2) = fun(x(2));
max1=1000; %最大迭代次数
for k = 2: max1
x(k+1) = x(k)-fun(x(k))*(x(k)-x(1))/(fun(x(k))-fun(x(1))); %两点弦割法的迭代公式
y(k+1) = fun(x(k+1));
e = abs(x(k+1)-x(k)); %误差区间
if e < err1
break
end
end
xx=x(k+1) %牛顿迭代法求得的近似解
figure;
plot([1:k+1],y)
xlabel('迭代步数');
ylabel('误差');
grid;
附录11 利用牛顿法求解非线性方程组
% 附录11 利用牛顿法求解非线性方程组
%% 利用牛顿法求解非线性方程组
clear %清除变量
clc %清除命令行窗口代码
format long
err1 = input('\n请输入最大允许误差：\n');
x0 = input('\n请输入给定初始向量(元素之间逗号间隔)：\n');
[allx,ally,xx,n]=NNewton(fun,x0,err1)
xx1=double(xx)
fid=fopen('data_xx.txt','w');
fprintf(fid,'%.8f\n',xx1'); %将方程组数值解存于'data_xx.txt'文件夹中
fclose(fid);
function [allx,ally,r,n]=NNewton(F,x0,eps)
if nargin==2 %函数调用中给定函数输入参数的数目为2
eps=1.0e-8;
end
x0 = transpose(x0); %transpose转置向量或矩阵
Fx = subs(F,transpose(symvar(F)),x0);
var = transpose(symvar(F));
dF = jacobian(F,var);
dFx = subs(dF,transpose(symvar(F)),x0);
n=dFx;
r=x0-inv(dFx)*Fx';
n=1;tol=1;
N=100;
symx=length(x0);
ally=zeros(symx,N);
allx=zeros(symx,N);
while tol>eps
x0=r;
Fx = subs(F,transpose(symvar(F)),x0);
dFx = subs(dF,transpose(symvar(F)),x0);
r=vpa(x0-inv(dFx)*Fx'); %牛顿法迭代
tol=norm(r-x0)
if(n>N)
disp('迭代步数太多，可能不收敛。');
break;
end
allx(:,n)=x0;
ally(:,n)=Fx;
n=n+1;
end
end
function f = fun(x) %非线性方程组，根据题目或实际情况进行更改
k=3; % k为所求解方程个数，同时也是自变量的个数
for i=1:k %将自变量变成符号
x(i)=sym (['x',num2str(i)]);
end
f(1)=x(1)^2+x(2)^2+x(3)^2-1;
f(2)=2*x(1)^2+x(2)^2-4*x(3)^2;
f(3)=3*x(1)^2-4*x(2)^2+x(3)^2;
end附录12 利用弦割法求解非线性方程组
% 附录12 利用弦割法求解非线性方程组
%% 利用弦割法求解非线性方程组(对应7.3第一个)
clear %清除变量
clc %清除命令行窗口代码
err1 = input('\n请输入最大允许误差：\n'); % 求解精度
x0 = input('\n请输入给定初始向量(元素之间逗号或空格间隔)：\n');
x0=x0';
format long
f = @(x)[x(1).^2+x(2).^2+x(3).^2-1; %待求函数，可以更改
2*x(1).^2+x(2).^2-4*x(3);
3*x(1).^2-4*x(2).^2+x(3).^2];
h = 0.01; % 步长
n=1; %迭代次数
z = [f(x0 + h*[1; 0; 0]), f(x0 + h*[0; 1; 0]), f(x0 + h*[0; 0; 1])]\f(x0);
x1 = x0 + h*z/(sum(z) - 1);
xx(:,n)=x1;
while norm(x1 - x0,2) > err1 %弦割法迭代求解
x0 = x1; n=n+1;
z = [f(x0 + h*[1; 0; 0]), f(x0 + h*[0; 1; 0]), f(x0 + h*[0; 0; 1])]\f(x0);
x1 = x0 + h*z/(sum(z) - 1);
xx(:,n)=x1;
end
csvwrite('data_xx.txt', x1);
%% 利用弦割法求解非线性方程组(对应7.3第二个)
clear %清除变量
clc %清除命令行窗口代码
err1 = input('\n请输入最大允许误差：\n'); % 求解精度
x0 = input('\n请输入给定初始向量(元素之间逗号或空格间隔)：\n');
x0=x0';
format long
f = @(x)[cos(x(1).^2+0.4*x(2))+x(1)^2+x(2)^2-1.6; %待求函数，可以更改
1.5*x(1).^2-1/0.36*x(2).^2-1.0];
h = 0.01; % 步长
n=1; %迭代次数
z = [f(x0 + h*[1; 0]), f(x0 + h*[0; 1]) ]\f(x0);
x1 = x0 + h*z/(sum(z) - 1);
xx(:,n)=x1;
while norm(x1-x0,2) > err1 %弦割法迭代求解
x0 = x1; n = n+1;
z = [f(x0 + h*[1; 0 ]), f(x0 + h*[0; 1])]\f(x0);
x1 = x0 + h*z/(sum(z) - 1);
xx(:,n)=x1;
end
csvwrite('data_xx.txt', x1);
附录13 利用布洛依登法求解非线性方程组
% 附录13 利用布洛依登法求解非线性方程组
%% 布洛依登法求解非线性方程组
% 注释部分对应第一组方程，非注释代码对应第二组方程
clear %清除变量
clc %清除命令行窗口代码
format long
err1 = input('\n请输入最大允许误差：\n');
x0 = input('\n请输入给定初始向量(元素之间逗号或空格间隔)：\n');
x0=x0';
[x,n,data]=broyden(x0,err1);
disp('计算结果为')
x
disp('迭代次数为')
n
data(:,1)=x0;
%抽取 data 中第一个变量数据画出曲线
subplot(2,1,1)
plot(data(1,:)),
title('x1 在迭代中的变化')
xlim([1.0 5.0]);ylim([0.90 1.10]);
%抽取 data 中的第二个变量数据 画出其变化曲线
subplot(2,1,2)
plot(data(2,:)),title('x2 在迭代中的变化')
xlim([1.0 5.0]);ylim([0.40 0.50]);
% subplot(3,1,1)
% plot(data(1,:)),
% title('x1 在迭代中的变化')
% xlim([1.0 16.0]);ylim([0.00 1.00]);
% %抽取 data 中的第二个变量数据 画出其变化曲线
% subplot(3,1,2)
% plot(data(2,:)),title('x2 在迭代中的变化')
% xlim([1.0 16.0]);ylim([0.00 1.00]);
%抽取第三个变量数据
% subplot(3,1,3)
% plot(data(3,:)),title('x3 在迭代中的变化')
% xlim([1.0 16.0]);ylim([0.00 1.00]);
%% 布洛依登法计算非线性方程组
function [x,n,data]=broyden(x0,tol)
%输入 x0 为迭代初值，tol 为误差容限 如果缺省 默认为 10 的-10 次方
%data 用来存放计算的中间数据便于计算收敛情况分析
if nargin==1
tol=1e-8;
end
% H0=df2(x0);
H0=df3(x0);
H0=inv(H0);
% x1=x0-H0*f2(x0);
x1=x0-H0*f3(x0);
n=1;
%设置初始误差 使之可以进入循环
wucha=0.1;
%循环迭代
while (wucha>tol)&(n<20) &(n<500)
wucha=norm(x1-x0);
dx=x1-x0;
% y=f2(x1)-f2(x0);
y=f3(x1)-f3(x0);
fenzi=dx'*H0*y; %是标量
H1=H0+(dx-H0*y)*(dx)'*H0/fenzi;
temp_x0=x0;
x0=x1;
% x1=temp_x0-H1*f2(temp_x0); %x1 的更新
x1=temp_x0-H1*f3(temp_x0); %x1 的更新
%更新 H 矩阵
H=H1;
n=n+1;
%data 用来存放中间数据
data(:,n)=x1;
end
x=x1;
end
% 目标方程组(对应7.3第一个)，注意方程组输出为列向量函数
function F = f2(x0)
x=x0(1);
y=x0(2);
z=x0(3);
ff1=x^2+y^2+z^2-1.0; %待求方程组，可更改
ff2=2*x^2+y^2-4*z;
ff3=3*x^2-4*y^2+z^2;
F=[ff1;ff2;ff3];
end
%Jacobi 矩阵,该方程输出为 n X n 维方阵(对应7.3第一个非线性方程组)
function f=df2(x0)
x=x0(1);
y=x0(2);
z=x0(3);
f=[2*x 2*y 2*z %目标方程组Jacobi矩阵，可以更改（第一列对x求导，第二列对y求导，依次类推）
4*x 2*y -4
6*x -8*y 2*z];
end
% 目标方程组(对应7.3第二个)，注意方程组输出为列向量函数
function F = f2(x0)
x=x0(1);
y=x0(2);
ff1=cos(x^2+0.4*y)+x^2+y^2-1.6; %待求方程组，可更改
ff2=1.5*x^2-1/0.36*y^2-1;
F=[ff1;ff2];
end
% Jacobi 矩阵,该方程输出为 2 X 2 维方阵(对应7.3第二个非线性方程组)
function f=df2(x0)
x=x0(1);
y=x0(2);
f=[-sin(x^2+0.4*y)*2*x+2*x -sin(x^2+0.4*y)*0.4+2*y %目标方程组Jacobi矩阵，可以更改（第一列对x求导，第二列对y求导）
3*x -2/0.36*y];
end
此处仅展示了代码实现，具体算法原理可参考《数值分析》等有关书籍，或在博主“资源”下载页面，下载相关文档进行查看。
扫一扫
4370
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
