【计算方法】部分习题的matlab求解

最新推荐文章于 2021-03-28 12:58:27 发布

加油加油再加油x

最新推荐文章于 2021-03-28 12:58:27 发布

阅读量3.3k

点赞数 18

分类专栏：数学建模文章标签： matlab

本文链接：https://blog.csdn.net/weixin_43325134/article/details/106354059

版权

数学建模专栏收录该内容

4 篇文章 3 订阅

订阅专栏

前言

要考试了整理一下。（没错，由于疫情原因，我们计算方法考试线下进行，可以用matlab。）

第一章、拉格朗日插值

在这里插入图片描述

clc,clear;
syms x;
% 拉格朗日插值点  x和y
x0 = [3,5,6];
y0 = [1.0986,1.6094,1.7918];
xq = 4;  %所求点的横坐标
len = length(x0);
%lx为核函数列表
lx = [];
for i = 1:len
	tmp = x0(:);
	tmp(i) = [];
	lxi = prod(x-tmp)/prod(x0(i)-tmp);
	lx = [lx lxi];
end
% y(x)为插值函数
y(x) = sum(lx.*y0);
disp("y = ");
disp(vpa(expand(y),4));
% 保留小数,求出结果
res = vpa(y(xq),4);
disp("y(xq)结果为")
disp(res);

%求误差的情况通常有一个已知的f（x）
f(x) = log(x);
% fx的n+1阶导数绝对值
fe(x) = abs(diff(f,len));

%尝试x，求最大值
x1 = x0(1):0.01:x0(len);
%求出n+1阶导数绝对值的最大值
fe0 = max(fe(x1));
%求出R(x)
R(x) = abs(fe0/factorial(len)*prod(x-x0));
disp("R(x)=")
disp(vpa(max(R(x1)),4))
%画出R（x）图像
% plot(x1,R(x1));

第二章、牛顿插值

在这里插入图片描述

clc,clear;
syms x;
%牛顿插值已知点
x0 = [3,4,5,6];
y0 = [1.0986,1.3863,1.6094,1.7918];
%所求值
xq = 4.5;
%求出x0的长度
len = length(x0);
%构建差商表table
table = zeros(len,len);
table(:,1) = y0';
for i=2:len
	tmp1 = table(i-1:len-1,i-1);
	tmp2 = table(i:len,i-1);
	xtmp = (x0(1:len-i+1)-x0(i:len))';
	table(i:len,i) = (tmp1-tmp2)./xtmp;
end
% 打印差商表
disp("差商表为")
disp(table)
xishu = diag(table);
N(x) = x;
for i = 1:len
	if (i==1)
		N(x) = N(x) + xishu(i);
	else
		N(x) = N(x) + xishu(i)*prod(x-x0(1:i-1));
	end
end
N(x) = N(x) - x;
%化简输出结果
vpa(simplify(N(x)),4)
vpa(N(xq),4)

%误差求法同拉格朗日
f(x) = log(x);
% fx的n+1阶导数
fe(x) = abs(diff(f,len));

%尝试x，求最大值
x1 = x0(1):0.01:x0(len);
%求出n+1阶导数绝对值的最大值
fe0 = max(fe(x1));
%求出R(x)
R(x) = abs(fe0/factorial(len)*prod(x-x0));
disp("xq loss")
disp(vpa(R(xq),4))
disp("y loss")
disp(vpa(max(R(x1)),4))
%画出R（x）图像
% plot(x1,R(x1));

第三章、线性插值&三次样条插值

线性插值

在这里插入图片描述

clc,clear;
syms x;
%分段线性插值
f(x) = 1/(1+x^2);
a = 2;
b = 5;
wucha = 0.01;

%求处达到精度的划分
% m为区间个数
% h为区间长度
ddf(x) = abs(diff(f,2));
x1 = a:0.01:b;
m2 = max(ddf(x1));
h = (8*wucha/m2)^0.5;
m = ceil((b-a)/h);
h = (b-a)/m;

%求出各个区段表达式
x0 = a:h:b;
y0 = f(x0);
y = [];
for i=1:m
	y = [y vpa((x-x0(i+1))/(x0(i)-x0(i+1))*y0(i) + (x-x0(i))/(x0(i+1)-x0(i))*y0(i+1),4)];
end
vpa(y,4)

三次样条插值

在这里插入图片描述

clc,clear;
%三次样条插值
%求达到误差精度的最大区间数量

syms x;
f(x) = log(x);
a = 3;
b = 6;
wucha = 0.0001;

df(x) = abs(diff(f,4));
m = max(df(a:0.01:b));

h = (wucha*384/5/m)^(0.25);
n = ceil((b-a)/h);
h = (b-a)/n;

n
h

第四章、最佳一致逼近

在这里插入图片描述

clc,clear;
%最佳一致逼近
syms x;
a = -1;
b = 1;
f(x) = exp(x);
cishu = 3; %逼近次数
k = 1:cishu+1;
x0 = (b-a)/2*cos((2.*k-1).*pi./(2.*length(k)))+(a+b)/2;
y0 = f(x0);

len = length(x0);
%lx为核函数列表
lx = [];
for i = 1:len
	tmp = x0(:);
	tmp(i) = [];
	lxi = prod(x-tmp)/prod(x0(i)-tmp);
	lx = [lx lxi];
end
% y(x)为插值函数
y(x) = sum(lx.*y0);
% 保留小数,求出结果
res = vpa(expand(y),4);
disp(res);

%计算误差
df(x) = abs(diff(f,cishu+1));
mdf = max(df(a:0.01:b));
wucha = 1/2^(cishu)/factorial(cishu+1)*mdf;
vpa(wucha,4)

第五章、最佳平方逼近和曲线拟合的最小二乘法

最佳平方逼近

在这里插入图片描述

clc,clear;
%最佳平方逼近
syms x;
%roux
rou(x) = x/x;
fai = [1 x x^2];
a = 0;
b = 1;
f(x) = exp(x);

[faix,faiy] = meshgrid(fai,fai);
%求内积
left = int(rou*faix.*faiy,a,b)
right = int(f.*fai,a,b)'
%矩阵求逆求出系数
xishu = inv(left)*right;

% vpa(xishu,4)
S(x) = sum(fai'.*xishu);
disp("polynomial is")
vpa(S(x),4)
%误差处理
%平方误差
pingfang = abs(int(f^2*rou,a,b)-int(S*f*rou,a,b));
disp("sqrt error")
vpa(pingfang,4)
%最大值误差
disp("max error")
zuidazhi(x) = abs(f-S);
vpa(max(zuidazhi(a:0.01:b)),4)

曲线拟合的最小二乘法

在这里插入图片描述

clc,clear;
%曲线拟合
syms x;
x0 = [-2 -1 0 1 2];
y0 = [-8.2 -1.2 0 1.25 8.05];
cishu = 2;

len = length(x0);
left = zeros(cishu+1,cishu+1);
xs = [];
for i = 1:cishu+1
	xs = [xs;x^(i-1)];
end

for i=1:cishu+1
	for j=1:cishu+1
		f(x) = x^(i+j-2);
		left(i,j) = sum(f(x0));
	end
end

right = zeros(cishu+1,1);
for i = 1:cishu+1
	right(i) = sum(y0.*x0.^(i-1));
end
%求出系数
a = inv(left)*right;
%表达式
f(x) = sum(a.*xs);
disp("polynomial is")
vpa(f(x),4)
%求出误差
e = sum((f(x0)-y0).^2)^0.5;
disp("error is")
vpa(e,4)

机械求积方法

在这里插入图片描述

clc,clear;
syms x;
% 构造积分精度
jingdu = 3;
% 原先的被积函数
f(x) = exp(x);
% 积分区间
a = -1;
b = 2;

func = [];
x0 = linspace(a,b,jingdu+1);
for i=1:jingdu+1
	func = [func x^(i-1)];
end

left = zeros(jingdu+1);
for i=1:jingdu+1
	f_i(x) = func(i);
	left(i,:) = f_i(x0);
end

right = int(func',a,b);

A = inv(left)*right;
disp("A=")
vpa(A,4)

%求得的积分结果
intres = sum(A'.*f(x0));
disp("求得结果为")
vpa(intres,4)
%求误差
%真实误差
acc = int(f,a,b);
acce = abs(acc - intres);
disp("真实误差为");
vpa(acce,4)
%误差估计
df(x) = diff(f,jingdu+1);
maxf = max(abs(df(a:0.01:b)));
e = maxf/factorial(jingdu+1)*int(abs(prod(x-x0)),a,b);
disp("误差估计为");
vpa(e,4)

复合求积公式

在这里插入图片描述

clc,clear;
%课堂作业1
%复化求积公式
syms x;
a = 0;
b = 4;
n= 4;
h=(b-a)/n;
f(x)=exp(x);
%求精确解
disp("精确解为")
acc=vpa(int(f,a,b));
disp(acc)
%复化梯形公式
i=a+h:h:b-h;
%复化梯形公式算公式如下
Tn=h/2*(f(a)+2*sum(f(i))+f(b));
disp("复化梯形公式的计算结果为")
disp(vpa(Tn,4))
%复化辛普森公式
i=a+h:h:b-h;
j=a+h/2:h:b-h/2;
%复化辛普森公式计算公式如下
Sn=h/6*(f(a)+4*sum(f(j))+2*sum(f(i))+f(b));
disp("复化辛普生公式的结算结果为")
disp(vpa(Sn,4))
%复化梯形公式的误差
% disp("复化梯形公式的实际误差为")
% acce1=abs(acc-vpa(Tn));
% disp(acce1)
disp("复化梯形公式的误差")
df(x) = diff(f,2);
Rt = abs((b-a)/12*h^2)*max(abs(df(a:0.01:b)));
disp(vpa(Rt,4))
%复化辛普生公式的误差
% disp("复化辛普生公式的实际误差为")
% acce2=abs(acc-vpa(Sn));
% disp(acce2)
disp("复化辛普生公式的误差为")
df(x) = diff(f,4);
Rs = abs((b-a)/2880*h^4)*max(abs(df(a:0.01:b)));
disp(vpa(Rs,4))

线性方程组的直接解法

在这里插入图片描述

clc,clear;
syms x;
%第一问不写了自个拿着计算器解去吧！
%判读是否收敛:求cond(A)
A = [8 -2 1
	 3 1 1.001
	 1 1 1];
% Ainv = inv(A);
%直接求出结果
cond(A,inf)
%经过一定过程
Ainv = inv(A);
A = abs(A);
Ainv =abs(Ainv);
max(sum(A,2))*max(sum(Ainv,2))
%当输出结果远远大于1时为病态的。

解线性方程组的迭代法

雅克比迭代法

在这里插入图片描述

clc,clear;
% 雅克比迭代
% 求解方法AX = B
syms x;
A = [4 1 2
	 1 4 1
	 2 -1 4];
B = [7;6;5];
len = length(B);
% 初值为0
X = zeros(len,1);
% 迭代次数
n = 3;
diagA = diag(A);
left = A;
left(logical(eye(size(left)))) = 0;
left = -left./diagA;
left0 = B./diagA;
for i = 1:n
	fprintf("第%d次迭代",i);
	%若不加分号显示迭代过程
	X = left * X + left0
end
% X
%B其实就是left
left

高斯-塞德尔迭代

在这里插入图片描述

clc,clear;
% 高斯赛德尔迭代
% 求解方法AX = B
syms x;
A = [4 1 2
	 1 4 1
	 2 -1 4];
B = [7;6;5];
len = length(B);
% 初值为0
X = zeros(len,1);
% 迭代次数
n = 3;

D = diag(diag(A));
L = -tril(A,-1);
U = -triu(A,1);
G = inv(D-L)*U;
d = inv(D-L)*B;
for i = 1:n
	fprintf("第%d次迭代",i);
	X = G * X + d
end
G

解线性方程组的迭代法的收敛性

雅克比迭代法

在这里插入图片描述

clc,clear;
%解线性方程组的收敛性
% 雅克比迭代（和老师答案有点不太一致）
% 求解方法AX = B
syms x;
A = [3 0 1
	 1 3 2
	 0 1 3];
B = [4;3;3];
jingdu = 0.001;

len = length(B);
diagA = diag(A);
left = A;
left(logical(eye(size(left)))) = 0;
left = -left./diagA;
left0 = B./diagA;
if(max(abs(eig(left)))>=1)
	disp("发散无法求出迭代次数");
end
X0 = zeros(len,1);
X1 = left * X0 + left0;

%left的范数
%我们主要看两个范数
%1.无穷范数
%2.F范数
condB = max(sum(abs(left),2)); %X0与X1的无穷范数
condX = max(abs(X1 - X0)); %无穷范数
disp("无穷范数求得迭代次数为")
k = ceil(log(jingdu/condX*(1-condB))/log(condB))
condB = sum(sum(left.^2))^0.5; %X0与X1的F范数
condX = sum((X1-X0).^2)^0.5; %F范数
disp("F范数求得迭代次数为")
k =  ceil(log(jingdu/condX*(1-condB))/log(condB))

高斯-塞德尔迭代

在这里插入图片描述

clc,clear;
% 解线性方程组的收敛性
% 高斯赛德尔迭代
% 求解方法AX = B
% 和老师的结果差1
syms x;
A = [3 0 1
	 1 3 2
	 0 1 3];
B = [4;3;3];
jingdu = 0.001;
len = length(B);
D = diag(diag(A));
L = -tril(A,-1);
U = -triu(A,1);
G = inv(D-L)*U;
d = inv(D-L)*B;

if(max(abs(eig(G)))>=1)
	disp("发散无法求出迭代次数");
end
X0 = zeros(len,1);
X1 = G * X0 + d;
%G的范数
%我们主要看两个范数
%1.无穷范数
%2.F范数
condB = max(sum(abs(G),2)); %X0与X1的无穷范数
condX = max(abs(X1 - X0)); %无穷范数
disp("无穷范数求得迭代次数为")
k = ceil(log(jingdu/condX*(1-condB))/log(condB))
condB = sum(sum(G.^2))^0.5; %X0与X1的F范数
condX = sum((X1-X0).^2)^0.5; %F范数
disp("F范数求得迭代次数为")
k =  ceil(log(jingdu/condX*(1-condB))/log(condB))

非线性方程的数值解法

二分法

在这里插入图片描述

clc,clear;
%非线性方程的数值解法
%二分法
%有可能还需要求在这个区间上有几个解，验证解的唯一性
syms x;
% 方程为f(x) = 0
f(x) = x^3 - x^2 + x - 1;
%迭代区间
a = 0;
b = 2;
%n为迭代次数
n = 5;
x1 = a;
x2 = b;
for i=1:n
	if (f((x1+x2)/2)*f(x1)>0)
		x1 = (x1+x2)/2;
	elseif(f((x1+x2)/2)*f(x2)>0)
		x2 = (x1+x2)/2;
	else
		x1 = (x1+x2)/2;
		x2 = (x1+x2)/2;
	end
end
disp("解在下列两个数之间")
x1,x2

不动点迭代法

在这里插入图片描述

clc,clear;
%非线性方程的数值解法
%不动点迭代法
syms x;
% 方程为f(x)为构造的不动点迭代方程
f(x) = (x^2 - x + 1)^(1/3);
%迭代区间
a = 0.5;
b = 1.5;
%n为迭代次数
n = 5;
%x0 为开始迭代的点
x0 = 0.5;

df(x) = diff(f,1);
% vpa(sum(abs(df(a:0.01:b))>1),5)
if(sum(logical(abs(df(a:0.01:b))>1))==0)
	disp("ok")
end
X = x0;
for i=1:n
	X = f(X);
end
vpa(X,4)

牛顿迭代法

在这里插入图片描述

clc,clear;
%非线性方程的数值解法
%牛顿迭代法
syms x;
% 方程为f(x) = 0
f(x) = x^3 - x^2 + x - 1;
%n为迭代次数
n = 5;
%x0 为开始迭代的点
x0 = 0.5;

df(x) = diff(f,1);
X = x0;
for i=1:n
	X = X - f(X)/df(X);
end
vpa(X,4)

常微分方程初值问题数值解法

欧拉法

clc,clear;
%欧拉法
h=0.05; %步长  
xf=2; %区间最右边x坐标
x = 0:h:2;%初始化x向量
y = zeros(1,length(x));%初始化y  
y(1) = 1; %y(0)=1  
F_xy = @(x, y)(y.^2-2*x)./y%y导数
for i=1:(length(x)-1) %循环迭代  
    y(i+1)=y(i) + h*(y(i)-2*x(i)/y(i))  
end  
% 	x  
% 	y  
%调用Matlab方法ode45()进行常微分方程初值问题的精确解求解  
%求得的解在函数调用的左端[x_ode45, y_ode45]  
x0 = x(1);  
y0 = y(1);  
yspan = [x0 xf];  
[x_ode45, y_ode45] = ode45(F_xy,yspan,y0);  
%画出精确解  
plot(x_ode45,y_ode45,'--');   
xlabel('x');ylabel('y');  
hold on;
%画出近似解  
plot(x,y,'-');
disp("欧拉法解为")
for i=1:length(y)
    fprintf("%.8f\n",y(i))
end
disp("精确解为")
for i=1:length(y_ode45)
    fprintf("%.8f\n",y_ode45(i));
end
xlabel('x');ylabel('y');  
legend('Exact','Euler');

改进欧拉法

clc,clear;
%改进欧拉法
h=0.05; %步长  
xf=2; %区间最右侧点  
x = 0:h:2;%初始化x向量  
y = zeros(1,length(x));%初始化y  
y(1) = 1; %y(0)=1  
F_xy = @(x, y)(y.^2-2*x)./y%y导数 
for i=1:(length(x)-1) %循环迭代  
    yp=y(i)+h*(y(i)-2*x(i)/y(i))  
    yc=y(i)+h*(yp-2*x(i+1)/yp)  
    y(i+1)=(yp+yc)/2  
end  
% 	x  
% 	y  
%调用Matlab方法ode45()进行常微分方程初值问题的精确解求解  
%求得的解在函数调用的左端[x_ode45, y_ode45]  
x0 = x(1);  
y0 = y(1);  
yspan = [x0 xf];  
[x_ode45, y_ode45] = ode45(F_xy,yspan,y0);  
%画出精确解  
plot(x_ode45,y_ode45,'--');   
xlabel('x');ylabel('y');
hold on;
%画出近似解  
plot(x,y,'-');  
y=vpa(y,8);
y_ode45=vpa(y_ode45,8);
    disp("改进欧拉法解为")
for i=1:length(y)
    fprintf("%.8f\n",y(i))
end
disp("精确解为")
for i=1:length(y_ode45)
    fprintf("%.8f\n",y_ode45(i));
end
xlabel('x');ylabel('y');  
legend('Exact','improve_Euler');

PS：csdn是不是没法插入matlab代码？我用的php的格式。
还有那个我求函数最大值是通过取很多点来实现的（为了省事），需要注意一下求出来的不一定是最大值。
以后可能会更新。
链接：https://pan.baidu.com/s/1W7YiEsf0KamzTyLxjTC-zg
提取码：blj5
可能有错误，各位大佬发现别忘记在评论区评论。
要是觉得有用，别忘记分享点赞哦！