下载文档原格式
数值分析中求解线性方程组的MATLAB程序（6种）1.回溯法（系数矩阵为上三角）function X=uptrbk(A,B)%求解方程组，首先化为上三角，再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A was singular.No unique solution.'break;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));2.系数矩阵为下三角function x=matrix_down(A,b)%求解系数矩阵是下三角的方程组n=length(b);x=zeros(n,1);x(1)=b(1)/A(1,1);for k=2:1:nx(k)=(b(k)-A(k,1:k-1)*x(1:k-1))/A(k,k);end3.普通系数矩阵（先化为上三角，在用回溯法）function X=uptrbk(A,B)%求解方程组，首先化为上三角，再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A was singular.No unique solution.'break;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));4.三角分解法function [X,L,U]=LU_matrix(A,B)%A是非奇异矩阵%AX=B化为LUX=B，L为下三角，U为上三角%程序中并没有真正解出L和U，全部存放在A中[N,N]=size(A);X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);R=1:N;for p=1:N-1[max1,j]=max(abs(A(p:N,p)));C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=R(p);R(p)=R(j+p-1);R(j+p-1)=d;if A(p,p)==0'A is singular.No unique solution'break;endfor k=p+1:Nmult=A(k,p)/A(p,p);A(k,p)=mult;A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);endendY(1)=B(R(1));for k=2:NY(k)=B(R(k))-A(k,1:k-1)*Y(1:k-1);endX(N)=Y(N)/A(N,N);for k=N-1:-1:1X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);endL=tril(A,-1)+eye(N)U=triu(A)5.雅克比迭代法function X=jacobi(A,B,P,delta,max1);%雅克比迭代求解方程组N=length(B);for k=1:max1for j=1:NX(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)breakendendX=X';k6.盖斯迭代法function X=gseid(A,B,P,delta,max1);%盖斯算法，求解赋初值的微分方程N=length(B);for k=1:max1for j=1:Nif j==1X(1)=(B(1)-A(1,2:N)*P(2:N))/A(1,1);elseif j==NX(N)=(B(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);elseX(j)=(B(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);endenderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)break;endendX=X';k。
数值分析（matlab程序）曹德欣曹璎珞第一章绪论数值稳定性程序，计算P11 试验题一积分function try_stableglobal nglobal aa=input('a=');N = 20;I0 = log(1+a)-log(a);I = zeros(N,1);I(1) = -a*I0+1;for k = 2:NI(k) = - a*I(k-1)+1/k;endII = zeros(N,1);if a>=N/(N+1)II(N) = (1+2*a)/(2*a*(a+1)*(N+1));elseII(N) =(1/((a+1)*(N+1))+1/N)/2;endfor k = N:-1:2II(k-1) = ( - II(k) +1/k) / a;endIII = zeros(N,1);for k = 1:Nn = k;III(k) = quadl(@f,0,1);endclcfprintf('\n 算法1结果算法2结果精确值') for k = 1:N,fprintf('\nI(%2.0f) %17.7f %17.7f %17.7f',k,I(k),II(k),III(k)) endfunction y = f(x)global nglobal ay = x.^n./(a+x);return第二章非线性方程求解下面均以方程y=x^4+2*x^2-x-3为例：1、二分法function y=erfen(a,b,esp)format longif nargin<3 esp=1.0e-4;endif fun(a)*fun(b)<0n=1;c=(a+b)/2;while c>espif fun(a)*fun(c)<0b=c;c=(a+b)/2;elseif fun(c)*fun(b)<0a=c;c=(a+b)/2;else y=c; esp=10000;endn=n+1;endy=c;elseif fun(a)==0y=a;elseif fun(b)==0y=b;else disp('these,nay not be a root in the intercal')endnfunction y=fun(x)y=x^4+2*x^2-x-3;2、牛顿法function y=newton(x0)x1=x0-fun(x0)/dfun(x0);n=1;while (abs(x1-x0)>=1.0e-4) & (n<=100000000)x0=x1;x1=x0-fun(x0)/dfun(x0);n=n+1;endy=x1nfunction y=fun(x)y=x^4+2*x^2-x-3; 3、割线法function y=gexian(x0,x1)x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %根据初始XO 和X1求X2 n=1;while (abs(x1-x0)>=1.0e-4) & (n<=100000000) %判断两个条件截止 x0=x1; %将x1赋给x0 x1=x2; %将x2赋给x1 x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %迭代运算 n=n+1; end y=x2 nfunction y=fun(x)y=x^4+2*x^2-x-3;第四章题目：推导外推样条公式：⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--------1232123211223322~~~~~22~n n n n n n n n d d d d M M M Mδμλμλμλδ，并编写程序与Matlab 的Spline 函数结果进行对比，最后调用追赶法解方程组。
2021数学建模C题代码详细MATLAB一、准备工作我们需要明确2021 数学建模C 题的具体要求和背景。
在这次比赛中，题目要求参赛选手利用给定的数据集和背景信息，通过建立数学模型，并编写相关代码进行模拟和求解。
针对不同的场景和要求，我们需要结合 MATLAB 软件的特点和优势，从数据处理、模型建立到结果分析，全面展现我们的解题思路和方法。
二、数据处理及可视化分析在开始编写 MATLAB 代码之前，我们需要对给定的数据集进行分析和处理。
具体来说，可以通过 MATLAB 的数据导入工具，将原始数据导入到 MATLAB 评台中，并进行数据预处理、清洗和可视化分析。
在处理数据的过程中，可以选用 MATLAB 中丰富的数据处理函数和绘图工具，展现数据的分布特点、趋势规律以及异常情况。
三、模型建立和求解在数据处理的基础上，我们将重点转向模型的建立和求解。
在C 题中，可能涉及到不同的数学模型，比如优化模型、拟合模型、动态模型等。
我们可以逐步引入相关的 MATLAB 函数和工具箱，比如优化工具箱、符号计算工具箱等，来构建模型并进行求解。
在代码的编写过程中，需要考虑到代码的优化和可读性，以便审阅人员和评审团队能够清晰地理解我们的解题思路和方法。
四、结果分析与讨论我们需要对模型的求解结果进行分析和讨论。
通过 MATLAB 的数据分析和统计工具，我们可以对求解出的参数和变量进行分析，验证模型的有效性和可行性。
我们也可以将结果以图表的形式进行展示，以便更直观地呈现我们的研究成果。
在结果分析的过程中，需要结合实际问题背景和应用场景，提出我们的见解和建议，从而全面、深刻地理解和解释模型的意义和作用。
五、个人观点和总结在这篇文章中，我个人的观点是，MATLAB 是一款功能强大的数学建模软件，具有丰富的数学函数和工具箱，能够满足不同场景和要求下的建模和求解需求。
通过编写 MATLAB 代码，不仅可以加深我们对数学模型和方法的理解，还能够培养我们的逻辑思维和问题解决能力。
1，秦九韶算法，求出P(x=3)=2+4x+5x^2+2x^3的值clear all;x=3;n=3;a(1)=2;a(2)=4;a(3)=5;a(4)=2v(1)=a(n+1);for k=2:(n+1);v(k)=x*v(k-1)+a(n-k+2);endp=v(n+1)p=,1132，一次线型插值程序：利用100.121.求115的开方。
clear all;x1=100;x2=121;y1=10;y2=11;x=115;l1=(x-x2)/(x1-x2);l2=(x-x1)/(x2-x1);p1=l1*y1+l2*y2p1=10.71433，分段插值程序，已知为S1（x)为（0，0），（1，1），（2，5）（3,8）上的分段一次插值，求S1（1.5）.clear allx=[0123];y=[0158];n=length(x);a=1.5;for i=2:nif(x(i-1)<=a<x(i));endendH1=y(i-1)+(y(i)-y(i-1))/(x(i)-x(i-1))*(a-x(i-1))H1=3.50004）曲线拟合：用一个5次多项式在区间[0，2π]内逼近函数sin(x)。
clear allX=linspace(0,2*pi,50);Y=sin(X);[P,S]=polyfit(X,Y,5)plot(X,Y,'k*',X,polyval(P,X),'k-')P=-0.00560.0874-0.39460.26850.87970.0102S=R:[6x6double]df:44normr:0.03375)求有理分式的导数clear allP=[3,5,0,-8,1,-5];Q=[10,5,0,0,6,0,0,7,-1,0,-100];[p,q]=polyder(P,Q)6）将以下数据按从小到大排序：4.3 5.7 5.2 1.89.4a=[4.35.75.21.89.4];b(1:100)=0;n=1;b(a*10)=1;for k=1:100a(n)=k/10;if b(k)>0a(n)=k/10;n=n+1;endendaa=1.8000 4.3000 5.2000 5.70009.400010.00007)用二分法求方程x 3-x-1=0在[1,2]内的近似根，要求误差不超过10-3。
一、拉格朗日插值#include<stdio.h>#include<stdlib.h>#include<math.h>void Lagrange(float s){double x[5]={0.2,0.4,0.6,0.8,1.0},y[5]={0.98,0.92,0.81,0.64,0.38},f,L=0;int i,j;for (i=0;i<5;i++){ f=1; for (j=0;j<5;j++) if(j!=i) f=(s-x[j])/(x[i]-x[j])*f; L+=f*y[i]; } printf("输出：%f\n",L);}void main(){float x;printf("输入插值点：");scanf("%f",&x);Lagrange(x);}二、牛顿插值#include<stdlib.h>#include<stdio.h>#include<math.h>int ND(float s){double x[5]={0.2,0.4,0.6,0.8,1.0},y[5]={0.98,0.92,0.81,0.64,0.38},p=0,g,f;int i,j,k;for (i=0;i<5;i++){for (j=4;j>i;j--){ f=x[j]-x[j-i-1];y[j]=(y[j]-y[j-1])/f;}g=y[i+1];for (k=0;k<=i;k++)g=g*(s-x[k]);p=p+g;}printf("输出插值点函数值：%f\n",p+y[0]);return 1;}void main(){float x;printf("输入插值点：");scanf("%f",&x);ND(x);}三、埃尔米特插值#include<stdio.h>#include<stdlib.h>#include<math.h>void Hermite(float s){double x[3]={0.25,1,2.25},y[3]={0.125,1,3.375},z[3]={0.75,1.5,2.25};double H=0.0,a,b,f,g;int i,j;for (i=0;i<3;i++){f=1.0;g=0.0;for (j=0;j<3&&j!=i;j++){f=f*(s-x[j])/(x[i]-x[j]);g=g+1/(x[i]-x[j]);}a=(1-2*(s-x[i])*g)*f*f;b=(s-x[i])*f*f;H=H+y[i]*a+z[i]*b;}printf("%f\n",H);}void main(){float x;printf("输入插值点：");scanf("%f",&x);Hermite(x);}四、三次样条插值#include <math.h>#include <stdio.h>#include <stdlib.h>void main(){int N=7,R=2,i,k;double p1,p2,p3,p4;double x[8]={0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9};double y[8]={0.4794,0.6442,0.7833,0.8912,0.9636,0.9975,0.9917,0.9463};double P0=-0.4794,Pn=-0.9463,u[3]={0.6,0.8,1.2},s[3];double h[7],a[8],c[7], g[8],af[8],ba[7],m[8];for(k=0;k <N;k++)h[k]=x[k+1]-x[k];for(k=1;k <N;k++) a[k]=h[k]/(h[k]+h[k-1]);for(k=1;k <N;k++) c[k]=1-a[k];for(k=1;k <N;k++) g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]);c[0]=a[N]=1;g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;ba[0]=c[0]/2; g[0]=g[0]/2;for(i=1;i <N;i++){ af[i]=2-a[i]*ba[i-1]; g[i]=(g[i]-a[i]*g[i-1])/af[i]; ba[i]=c[i]/af[i]; }af[N]=2-a[N]*ba[N-1];g[N]=(g[N]-a[N]*g[N-1])/af[N];m[N]=g[N];for(i=N-1;i>=0;i--) m[i]=g[i]-ba[i]*m[i+1];for(i=0;i <=R;i++){ k=0;while(u[i]> x[k+1])k++;p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3);p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3);p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2);p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2);s[i]=p1+p2+p3+p4;}printf( "\nx= ");for(i=0;i <=N;i++)printf( "%8.1f ",x[i]);printf( "\ny= ");for(i=0;i <=N;i++)printf( "%8.4f ",y[i]);printf( "\n\nu= ");for(i=0;i <=R;i++)printf( "%9.2f ",u[i]);printf( "\n插值点：s= ");for(i=0;i <=R;i++)printf( "%9.5f ",s[i]);printf("\n");}五、复合梯形公式#include<stdio.h>#include<stdlib.h>#include<math.h>double FTX(int n,float a,float b){double f=0,t,h,*x,*y;int i;x=(double*)malloc((n+1)*sizeof(double));y=(double*)malloc((n+1)*sizeof(double));h=(b-a)/n;for(i=0;i<n+1;i++){x[i]=a+i*h;y[i]=sin(x[i]);}for(i=1;i<n;i++) f=f+2*y[i];t=h/2*(y[0]+f+y[n]);printf("输出函数值：%f\n",t);return 1;}void main(){float a,b;int n;printf("输入区间上，下限：");scanf("%f %f",&a,&b);printf("输入等分区间数：");scanf("%d",&n);FTX(n,a,b);}六、复合辛普森求积公式#include<stdio.h>#include<stdlib.h>#include<math.h>double FSP(int n,float a,float b){double f1=0,f2=0,h,*x1,*y1,*x2,*y2;int i;x1=(double*)malloc((n+1)*sizeof(double));y1=(double*)malloc((n+1)*sizeof(double));x2=(double*)malloc(n*sizeof(double));y2=(double*)malloc(n*sizeof(double));h=(b-a)/n;for(i=0;i<n+1;i++){x1[i]=a+h*i;y1[i]=sin(x1[i])*x1[i];}for(i=0;i<n;i++){x2[i]=x1[i]+h/2;y2[i]=sin(x2[i])*x2[i];}for(i=1;i<n;i++) f1=f1+2*y1[i];for(i=0;i<n;i++) f2=f2+4*y2[i];printf("输出函数值：%f\n",h/6*(y1[0]+f1+f2+y1[n]));return 1;}void main(){float a,b;int n;printf("输入区间上，下限：");scanf("%f %f",&a,&b);printf("输入等分区间数：");scanf("%d",&n);FSP(n,a,b);}七、直接三角分解法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){doubleA[3][3]={0.25,0.2,0.166667,0.3333,0.25,0.2,0.5,1,2},x[3],y[3],b[3]={9,8,8},L[3][3],U[3][3],f1=0,f2=0;int i,j,k;for(i=0;i<3;i++)for(j=0;j<3;j++){U[i][j]=0;L[i][j]=0;}for(i=0;i<3;i++){U[0][i]=A[0][i];L[i][0]=A[i][0]/U[0][0];L[i][i]=1;}for(i=1;i<3;i++)for(j=i;j<3;j++){for(k=0;k<=i-1;k++){f1=f1+L[i][k]*U[k][j];f2+=L[j][k]*U[k][i];} U[i][j]=A[i][j]-f1;L[j][i]=(A[j][i]-f2)/U[i][i];f1=0;f2=0;}y[0]=b[0];for(i=1;i<3;i++){for(j=0;j<=i-1;j++) f1+=L[i][j]*y[j];y[i]=b[i]-f1;f1=0;}x[2]=y[2]/U[2][2];for(i=1;i>=0;i--){for(j=i+1;j<3;j++) f2+=U[i][j]*x[j];x[i]=(y[i]-f2)/U[i][i];f2=0;} printf("输出L矩阵：\n");for(i=0;i<3;i++){for(j=0;j<3;j++)printf("%f ",L[i][j]);printf("\n");}printf("输出U矩阵：\n");for(i=0;i<3;i++){for(j=0;j<3;j++)printf("%f ",U[i][j]);printf("\n");}printf("输出求解结果：\n");for(i=0;i<3;i++)printf("%f ",x[i]);printf("\n");}八、改进的平方法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){double A[3][3]={2,-1,1,-1,-2,3,1,3,1},x[3],y[3],b[3]={4,5,6},d[3],L[3][3],U[3][3],f1=0,f2=0; int i,j,k,n=3;for(i=0;i<n;i++)for(j=0;j<n;j++) {U[i][j]=0;L[i][j]=0;}d[0]=A[0][0];L[0][0]=1;for(i=1;i<n;i++){L[i][i]=1;for(j=0;j<=i-1;j++) {for(k=0;k<=j-1;k++) f1+=U[i][k]*L[j][k];U[i][j]=A[i][j]-f1;L[i][j]=U[i][j]/d[j];f2+=U[i][j]*L[i][j];f1=0;}d[i]=A[i][j]-f2;f2=0;}y[0]=b[0];for(i=1;i<n;i++) {for(j=0;j<=i-1;j++) f1+=L[i][j]*y[j];y[i]=b[i]-f1;f1=0;}x[n-1]=y[n-1]/d[n-1];for(i=n-2;i>=0;i--) {for(j=i+1;j<n;j++) f2+=L[j][i]*x[j];x[i]=y[i]/d[i]-f2;f2=0;}printf("输出L矩阵：\n");for(i=0;i<n;i++){for(j=0;j<n;j++) printf("%f ",L[i][j]); printf("\n");}printf("输出U矩阵：\n");for(i=0;i<n;i++){for(j=0;j<n;j++) printf("%f ",U[i][j]); printf("\n");}printf("输出求解结果：\n");for(i=0;i<n;i++) printf("%f ",x[i]);printf("\n");}九、追赶法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){double A[5][5]={2,-1,0,0,0,-1,2,-1,0,0,0,-1,2,-1,0,0,0,-1,2,-1,0,0,0,-1,2},f[5]={1,0,0,0,0}; double x[5],y[5],U[5][5];int i,j,n=5;for(i=0;i<n;i++)for(j=0;j<n;j++) U[i][j]=0;U[0][0]=U[n-1][n-1]=1;U[0][1]=A[0][1]/A[0][0];for(i=1;i<n-1;i++){ U[i][i]=1;U[i][i+1]=A[i][i+1]/(A[i][i]-A[i][i-1]*U[i-1][i]);}y[0]=f[0]/A[0][0];for(i=1;i<n;i++)y[i]=(f[i]-A[i][i-1]*y[i-1])/(A[i][i]-A[i][i-1]*U[i-1][i]);x[n-1]=y[n-1];for(i=n-2;i>=0;i--) x[i]=y[i]-U[i][i+1]*x[i+1];printf("输出U矩阵：\n");for(i=0;i<n;i++){ for(j=0;j<n;j++) printf("%f ",U[i][j]); printf("\n");}printf("输出求解结果：\n");for(i=0;i<n;i++) printf("%f ",x[i]);printf("\n");}十、雅可比迭代法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[50][3],b[3]={-12,20,3},f=0,L=1;int n=3,i,j,k=1;for (i=0;i<3;i++) x[0][i]=0;/*初始值*/while(L<0.003){for(i=0;i<n;i++){ for(j=0;j<n;j++)if(j!=i) f=f+A[i][j]*x[k-1][j];x[k][i]=(b[i]-f)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i<n;i++)if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1) L=x[k][i]-x[k-1][i];k++;}printf("输出每次迭代求得的X值：\n");for(i=1;i<k;i++){printf("第%d次迭代：",i);for(j=0;j<n;j++)printf("%f ",x[i][j]);printf("\n");}printf("\n输出迭代次数：%d\n",k-1);}十一、高斯—塞德尔迭代法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[50][3],b[3]={-12,20,3},f1=0,f2=0,L=1;int i,j,k=1,n=3;for (i=0;i<3;i++) x[0][i]=0;while(L>0.00001||L<-0.00001){for(i=0;i<n;i++){ for(j=0;j<n;j++){if(j<i) f1+=A[i][j]*x[k][j];else if(j>i) f2+=A[i][j]*x[k-1][j];}x[k][i]=(b[i]-f1-f2)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i<n;i++)if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1) L=x[k][i]-x[k-1][i];k++;}printf("输出迭代X矩阵：\n");for(i=1;i<k;i++){for(j=0;j<n;j++) printf("%f ",x[i][j]); printf("\n");}printf("\n输出迭代次数：%d\n",k-1);}十二、超松弛迭代法#include<stdio.h>#include<stdlib.h>#include<math.h>void main(){double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[30][3],b[3]={-12,20,3},f1=0,f2=0,L=1,w=0.9;int i,j,k=1,n=3;for(i=0;i<3;i++) x[0][i]=0;while(L>0.0001||L<-0.0001){for(i=0;i<n;i++){ for(j=0;j<n;j++){if(j<i) f1+=A[i][j]*x[k][j];else if(j>i) f2+=A[i][j]*x[k-1][j];}x[k][i]=w*(b[i]-f1-f2)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i<n;i++)if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1) L=x[k][i]-x[k-1][i];k++;}printf("输出迭代X矩阵：\n");for(i=1;i<k;i++){for(j=0;j<n;j++) printf("%f ",x[i][j]); printf("\n");}printf("\n输出迭代次数：%d\n",k-1);}数值分析算法程序班级：计算08Q2班姓名：甄彦福。
列主元法function lianzhuyuan(A,b)n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵Ab=zeros(n,1); %矩阵bX=zeros(n,1); %解Xfor i=1:nfor j=1:nA(i,j)=(1/(i+j-1)); %生成hilbert矩阵Aendb(i,1)=sum(A(i,:)); %生成矩阵bendfor i=1:n-1j=i;top=max(abs(A(i:n,j))); %列主元k=j;while abs(A(k,j))~=top %列主元所在行k=k+1;endfor z=1:n %交换主元所在行a1=A(i,z);A(i,z)=A(k,z);A(k,z)=a1;enda2=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:nA(s,p)=A(s,p)-m*A(i,p);endb(s,1)=b(s,1)-m*b(i,1);endendX(n,1)=b(n,1)/A(n,n); %回代开始for i=n-1:-1:1s=0; %初始化sfor j=i+1:ns=s+A(i,j)*X(j,1);endX(i,1)=(b(i,1)-s)/A(i,i);endX欧拉法clcclear% 欧拉法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<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代过程for t=0:h:TY=(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;endt=[0:h:T];% 红线为解析解，'x'为数值解plot(t,S1,'r',t,S3,'x')改进欧拉法clcclearp=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<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代过程for t=0:h:TY=(exp(-t));R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t);k1=-y1;l1=y1-p*r1;k2=-(y1+h*k1);l2=y1+h*k1-p*(r1+h*l1);y=y1+h*(k1+k2)/2;r=r1+h*(k1+k2)/2;r1=r;y1=y;S1(i)=Y;S2(i)=R;S3(i)=y;S4(i)=r;i=i+1;endt=[0:h:T];% 红线为解析解，'x'为数值解plot(t,S1,'r',t,S3,'x')高斯-赛德尔function gaosisaideern=input('n='); %阶数tol=input('tol='); %迭代精度A=zeros(n,n);b=zeros(n,1); %生成b向量for i=1:n %给Hilbert矩阵和b向量赋值for j=1:nA(i,j)=(1/(i+j-1));endb(i,1)=sum(A(i,:));endy=zeros(n,1); %迭代解x1=zeros(n,1); %准确解for i=1:ny(i,1)=0; %迭代解赋初值x1(i,1)=1; %生成准确解endk=0;while norm(y-x1)>=tol %精度控制(采用自动步数控制) k=k+1;for i=1:n %迭代开始a1=0;a2=0;for j=1:i-1a1=a1+A(i,j)*y(j,1);endfor j=i+1:na2=a2+A(i,j)*y(j,1);endy(i,1)=(b(i,1)-a1-a2)/A(i,i);endenddisp('迭代步数k')kdisp(y) %显示yend最速下降法function gaosisaideern=input('阶数n='); %阶数tol=input('迭代精度tol='); %迭代精度eps=input('最速下降法eps=');A=zeros(n,n);b=zeros(n,1); %生成b向量for i=1:n %给Hilbert矩阵和b向量赋值for j=1:nA(i,j)=(1/(i+j-1));endb(i,1)=sum(A(i,:));endy=zeros(n,1); %迭代解x1=zeros(n,1); %准确解t=zeros(n,1);r=zeros(n,1);for i=1:ny(i,1)=0; %迭代解赋初值x1(i,1)=1; %生成准确解endr=b-A*y;while norm(r)>=eps; %先进行最速下降法求得进行赛德尔迭代的初始解yt=(r'*r)/(r'*A*r);s1=t*r;y=y+s1;r=b-A*y;endk=0;while norm(y-x1)>=tol %精度控制(采用自动步数控制)k=k+1;for i=1:n %迭代开始a1=0;a2=0;for j=1:i-1a1=a1+A(i,j)*y(j,1);endfor j=i+1:na2=a2+A(i,j)*y(j,1);endy(i,1)=(b(i,1)-a1-a2)/A(i,i);endenddisp('迭代步数k')disp(k)disp(y) %显示y四阶龙格-库塔法clcclearp=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<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代过程for t=0:h:TY=(exp(-t));R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t);k1=-y1;l1=y1-p*r1;k2=-(y1+h*k1/2);l2=y1+h*k1/2-p*(r1+h*l1/2);k3=-(y1+h*k2/2);l3=y1+h*k2/2-p*(r1+h*l2/2);k4=-(y1+h*k3);l4=y1+h*k3-p*(r1+h*l3);y=y1+h*(k1+2*k2+2*k3+k4)/6;r=r1+h*(l1+2*l2+2*l3+l4)/6;r1=r;y1=y;S1(i)=Y;S2(i)=R;S3(i)=y;S4(i)=r;i=i+1;endt=[0:h:T];% 红线为解析解，'x'为数值解plot(t,S1,'r',t,S3,'x')。
MATLAB 编程题库 1.下面的数据表近似地满足函数21cxbax y ++=，请适当变换成为线性最小二乘问题，编程求最好的系数c b a ,,，并在同一个图上画出所有数据和函数图像.625.0718.0801.0823.0802.0687.0606.0356.0995.0628.0544.0008.0213.0362.0586.0931.0ii 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.mfunction [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.mfunction v=f(x)v=x.*log(x)-1;>> edit g.mfunction z=g(y)z=y.^5+y-1;>> [x1 iter1]=gexianfa('f',1,3,1e-10)x1 =1.7632iter1 =6>> [x2 iter2]=gexianfa('g',0,1,1e-10)x2 =0.7549iter2 =83.使用GS 迭代求解下述线性代数方程组：123123123521242103103x x x x x x x x x ì++=-ïïïï-++=íïïï-+=ïî解：>> edit gsdiedai.mfunction [x iter]=gsdiedai(A,x0,b,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); iter=0; x=x0;while((norm(b-A*x)./norm(b))>tol) iter=iter+1; x0=x;x=(D-L)\(U*x0+b); end>> A=[5 2 1;-1 4 2;1 -3 10]; >> b=[-12 10 3]'; >>tol=1e-4; >>x0=[0 0 0]';>> [x iter]=gsdiedai(A,x0,b,tol); >>x x =-3.0910 1.2372 0.9802 >>iter iter = 64.用四阶Range-kutta 方法求解下述常微分方程初值问题（取步长h=0.01）,(1)2x dy y e xy dx y ìïï=++ïíïï=ïî解：>> edit ksf2.mfunction v=ksf2(x,y) v=y+exp(x)+x.*y;>> a=1;b=2;h=0.01; >> n=(b-a)./h; >> x=[1:0.01:2]; >>y(1)=2;>>fori=2:(n+1)k1=h*ksf2(x(i-1),y(i-1));k2=h*ksf2(x(i-1)+0.5*h,y(i-1)+0.5*k1); k3=h*ksf2(x(i-1)+0.5*h,y(i-1)+0.5*k2); k4=h*ksf2(x(i-1)+h,y(i-1)+k3); y(i)=y(i-1)+(k1+2*k2+2*k3+k4)./6; end >>y调用函数方法>> edit Rangekutta.mfunction [x y]=Rangekutta(f,a,b,h,y0) x=[a:h:b]; n=(b-a)/h; y(1)=y0; fori=2:(n+1)k1=h*(feval(f,x(i-1),y(i-1)));k2=h*(feval(f,x(i-1)+0.5*h,y(i-1)+0.5*k1)); k3=h*(feval(f,x(i-1)+0.5*h,y(i-1)+0.5*k2)); k4=h*(feval(f,x(i-1)+h,y(i-1)+k3)); y(i)=y(i-1)+(k1+2*k2+2*k3+k4)./6; end>> [x y]=Rangekutta('ksf2',1,2,0.01,2); >>y5.取0.2h =，请编写Matlab 程序，分别用欧拉方法、改进欧拉方法在12x ≤≤上求解初值问题。
1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6)disp('牛顿法');i=1;n0=180;p0=pi/3;tol=10^(-6);for i=1:n0p=p0-(p0-sin(p0))/(1-cos(p0));if abs(p-p0)<=10^(-6)disp('用牛顿法求得方程的根为')disp(p);disp('迭代次数为：')disp(i)break;endp0=p;endif i==n0&&~(abs(p-p0)<=10^(-6))disp(n0)disp('次牛顿迭代后无法求出方程的解')end2、disp('Steffensen加速');p0=pi/3;for i=1:n0p1=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;endp0=p;endif i==n0&&~(abs(p-p0)<=10^(-6))disp(n0)disp('次Steffensen加速后无法求出方程的解')end1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，%误差限为 e=10^-4disp('二分法')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:n0p=(a+b)/2;fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;if fp==0||(abs((b-a)/2)<tol)disp('用二分法求得方程的根p=')disp(p)disp('二分迭代次数为：')disp(i)break;endif fa*fp>0a=p;else b=p;endendif i==n0&&~(fp==0||(abs((b-a)/2)<tol))disp(n0)disp('次二分迭代后没有求出方程的根')end2、%使用牛顿法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，%误差限为 e=10^-4disp('牛顿法')p0=0.3;for i=1:n0p=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)<tol)disp('用牛顿法求得方程的根p=')disp(p)disp('牛顿迭代次数为：')disp(i)break;endp0=p;endif i==n0&&~(abs(p-p0)<tol)disp(n0)disp('次牛顿迭代后没有求出方程的根')end3、%使用割线法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，%误差限为 e=10^-4disp('割线法')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:n0p=p1-q1*(p1-p0)/(q1-q0);if abs(p-p1)<toldisp('用割线法求得方程的根p=')disp(p)disp('割线法迭代次数为：')disp(i)break;endp0=p1;q0=q1;pp=p1;p1=p;q1=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;endif i==n0&&~(abs(p-pp)<tol)disp(n0)disp('次割线法迭代后没有求出方程的根')end4、%使用试位法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，%误差限为 e=10^-4disp('试位法')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:n0p=p1-q1*(p1-p0)/(q1-q0);if abs(p-p1)<toldisp('用试位法求得方程的根p=')disp(p)disp('试位法迭代次数为：')disp(i)break;endq=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1;if q*q1<0p0=p1;q0=q1;endpp=p1;p1=p;q1=q;endif i==n0&&~(abs(p-pp)<tol)disp(n0)disp('次试位法迭代后没有求出方程的根')end5、%使用muller方法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根，%误差限为 e=10^-4disp('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:n0b=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)<abs(d+D))E=b+D;else E=b-D;endh=-2*(600*(x2.^4)-550*(x2.^3)+200*(x2.^2)-20*x2-1)/E;p=x2+h;if abs(h)<toldisp('用muller方法求得方程的根p=')disp(p)disp('muller方法迭代次数为：')disp(i)break;endx0=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);endif i==n0%条件有待商榷？！disp(n0)disp('次muller方法迭代后没有求出方程的根')end1、%观察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+1q(1,i)=y(i);endh=0.05;z=-1:h:1;for k=1:2/h+1for i=2:n+1for j=2:iq(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));endendw(k)=q(n+1,n+1);endhold 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+1q(1,i)=y(i);endh=0.05;z=-1:h:1;for k=1:2/h+1for i=2:n+1for j=2:iq(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));endendw(k)=q(n+1,n+1);endhold 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-1h(i)=x(i+1)-x(i);endA=[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:8b(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));endfor i=1:8z=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:6h(i)=x(i+1)-x(i);endA=[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:6b(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));endfor i=1:6z=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:3h(i)=x(i+1)-x(i);endA=[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:3b(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));endfor i=1:3z=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;endgrid on,title('注：横纵坐标的比例不一样！！！');1、%用不动点迭代法求方程 x-e^x+4=0的正根与负根，误差限是10^-6%disp('不动点迭代法');n0=100;p0=-5;for i=1:n0p=exp(p0)-4;if abs(p-p0)<=10^(-6)if p<0disp('|p-p0|=')disp(abs(p-p0))disp('不动点迭代法求得方程的负根为:')disp(p);break;elsedisp('不动点迭代法无法求出方程的负根.')endelsep0=p;endendif i==n0disp(n0)disp('次不动点迭代后无法求出方程的负根')endp1=1.7;for i=1:n0pp=exp(p1)-4;if abs(pp-p1)<=10^(-6)if pp>0disp('|p-p1|=')disp(abs(pp-p1))disp('用不动点迭代法求得方程的正根为')disp(pp);elsedisp('用不动点迭代法无法求出方程的正根');endbreak;elsep1=pp;endendif i==n0disp(n0)disp('次不动点迭代后无法求出方程的正根')end2、%用牛顿法求方程 x-e^x+4=0的正根与负根，误差限是10^-6 disp('牛顿法')n0=80;p0=1;for i=1:n0p=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;elsep0=p;endendif i==n0disp(n0)disp('次牛顿迭代后无法求出方程的解')endp1=-3;for i=1:n0p=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;elsep1=p;endendif i==n0disp(n0)disp('次牛顿迭代后无法求出方程的解')end1、使用欧拉法、改进欧拉法和四阶R-K方法求下列微分方程的解。
本文档包含上一个文档中的五个数值分析实验题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:7t=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.00001break;endend实验二C程序#include"stdio.h"#include"math.h"//已知量double x[10]={1,2,3,4,5,6,7,8,9,10};double fd1=1,fd2=0.1;double fx[10]={0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851}; //函数声明void canshu(double *h,double *r,double *u,double *d);void wanju(double *u,double *r,double *d,double *m);double yths(double *m,double *h,double p);double ytdhs(double *m,double *h,double p);//主函数void main(){double p,q1,q2;double h[9],u[9],r[9],m[10],d[10];// 矩阵Am=d参数数据u[8]=1,r[0]=1;canshu(h,r,u,d);wanju(u,r,d,m);printf("请输入待求数值：x=");scanf("%lf",&p);q1=yths(m,h,p);q2=ytdhs(m,h,p);printf("输出结果为:\nf(%lf)=%lf\nf'(%lf)=%lf\n",p,q1,p,q2);}//求参函数void canshu(double *h,double *r,double *u,double *d) {inti,j;double a[10];for(i=0;i<9;i++){h[i]=*(x+i+1)-*(x+i);a[i]=(fx[i]-fx[i+1])/(x[i]-x[i+1]);j=i-1;if(j>-1){r[i]=h[j+1]/(h[j+1]+h[j]);u[j]=h[j]/(h[j]+h[j+1]);d[i]=(6/(h[j]+h[i]))*(a[i]-a[j]);}}d[0]=(6/h[0])*(a[i-9]-fd1);d[9]=(6/h[i-1])*(fd2-a[i-1]);}//追赶法求弯矩向量mvoid wanju(double *u,double *r,double *d,double *m){inti,j;double l[9],n[10],z[10];n[0]=2;z[0]=d[0];for(i=0;i<9;i++){l[i]=u[i]/n[i];n[i+1]=2-l[i]*r[i];z[i+1]=d[i+1]-l[i]*z[i];}m[9]=z[9]/n[9];for(j=8;j>-1;j--){m[j]=(z[j]-r[j]*m[j+1])/n[j];}}//三次样条插值函数double yths(double *m,double *h,double p) {inti,j=0;double q;for(i=0;i<10;i++){if(x[i]<=(p-1)&&x[i+1]>=(p-1))break;}q=pow(x[i+1]-p,3)/6/h[i+1]*m[i]+pow(p-x[i],3)/6/h[i+1]*m[i+1]+(fx[i]-pow(h[i+1],2)*m[i]/6)*(x[i+1]-p)/h[i+1]+(fx[i+1]-pow(h[i+1],2)*m[i+1]/6)*(p-x[i])/h[i+1];return q;}//三次样条插值导函数double ytdhs(double *m,double *h,double p){inti,j=0;double q;for(i=0;i<10;i++){if(x[i]<=(p-1)&&x[i+1]>=(p-1))break;}q=-pow(x[i+1]-p,2)/2/h[i+1]*m[i]+pow(p-x[i],2)/2/h[i+1]*m[i+1]- (fx[i]-pow(h[i+1],2)*m[i]/6)/h[i+1]+(fx[i+1]-pow(h[i+1],2)*m[i+1]/6)/h[i+1];return q;}Matlab程序function[h,r,u,d]=canshu(x,fx,h,r,u,d,fd1,fd2)%UNTITLED10 Summary of this function goes here % Detailed explanation goes herefor i=1:9h(i)=x(i+1)-x(i);a(i)=(fx(i)-fx(i+1))/(x(i)-x(i+1));j=i-1;if j>0r(i)=h(j+1)/(h(j+1)+h(j));u(j)=h(j)/(h(j)+h(j+1));d(i)=(6/(h(j)+h(i)))*(a(i)-a(j));endendd(1)=(6/h(1))*(a(i-8)-fd1);d(10)=(6/h(i))*(fd2-a(i));endglobal x fd1 fd2 fx h d m r u;x=[1,2,3,4,5,6,7,8,9,10];fd1=1;fd2=0.1;fx=[0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.19 72246,2.3025851];u(9)=1;r(1)=1;[h,r,u,d]=canshu(x,fx,h,r,u,d,fd1,fd2);[m]=wanju(u,r,d,m);p=input('请输入待求数值:x=');q1=yths(m,h,p,x,fx);q2=ytdhs(m,h,p,fx,x);disp(['输出结果为:','f(',num2str(p),')=',num2str(q1),',','f(',num2str(p),')=',num2str(q2)]); function[m]=wanju(u,r,d,m)%UNTITLED12 Summary of this function goes here% Detailed explanation goes heren(1)=2;z(1)=d(1);for i=1:9l(i)=u(i)/n(i);n(i+1)=2-l(i)*r(i);z(i+1)=d(i+1)-l(i)*z(i);endm(10)=z(10)/n(10);for j=9:-1:1m(j)=(z(j)-r(j)*m(j+1))/n(j);endendfunction [q] =ytdhs(m,h,p,fx,x)%UNTITLED14 Summary of this function goes here% Detailed explanation goes herefor i=1:10if x(i)<=(p)&&x(i+1)>=(p)break;endq=-(x(i+1)-p)^2/2/h(i+1)*m(i)+(p-x(i))^2/2/h(i+1)*m(i+1)-(fx(i)-(h(i+1))^2*m(i)/6)/h(i+1)+(fx(i +1)-(h(i+1))^2*m(i+1)/6)/h(i+1);endendfunction[q]=yths(m,h,p,x,fx)%UNTITLED13 Summary of this function goes here% Detailed explanation goes herefor i=1:10if x(i)<=(p)&&x(i+1)>=(p)break;endq=(x(i+1)-p)^3/6/h(i+1)*m(i)+(p-x(i))^3/6/h(i+1)*m(i+1)+(fx(i)-(h(i+1))^2*m(i)/6)*(x(i+1)-p)/h(i+1)+(fx(i+1)-(h(i+1))^2*m(i+1)/6)*(p-x(i))/h(i+1); endend实验三C程序#include"stdio.h"#include"math.h"//已知量double b=3,a=1;//函数声明double romberg();double fhs(double x);//主函数void main(){ double q;q=romberg();printf("输出结果为：r=%lf\n",q);}//Romberg算法double romberg(){intn,j,k=0,m=0,z=0,i;double e=1,qh=0,q,h,x0,r1=0,t[50],s[50],c[50],r[50]; t[0]=(b-a)/2*(fhs(a)+fhs(b));for(i=1;1;i++){n=pow(2,i-1);h=(b-a)/n;x0=a+0.5*h;for(j=0;j<n;j++){qh+=fhs(x0+j*h);}t[i]=0.5*(t[i-1]+h*qh);qh=0;if(i>=3){for(k;k<i;k++)s[k]=(4*t[k+1]-t[k])/3;for(m;m<k-1;m++)c[m]=(16*s[m+1]-s[m])/15;for(z;z<m-1;z++)r[z]=(64*c[z+1]-c[z])/63;e=fabs(r1-r[z-1]);r1=r[z-1];}if(e<0.00001)break;}q=r[z-1];return q;}//f函数double fhs(double x){double y;y=pow(3,x)*pow(x,14)*(5*x+7)*sin(pow(x,2));return y;}Matlab程序function [y] =fhs(x)%UNTITLED4 Summary of this function goes here % Detailed explanation goes herey=3^x*x^14*(5*x+7)*sin(x^2);endfunction[q]=romberg()%UNTITLED3 Summary of this function goes here % Detailed explanation goes hereglobal b a n j k m z i e qh h x0;t=zeros(1,20);s=zeros(1,20);c=zeros(1,20);r=zeros(1,20);k=1;m=1;z=1;e=1;r1=0;t(1)=(b-a)/2*(fhs(a)+fhs(b));qh=0;for i=1:20n=2^(i-1);h=(b-a)/n;x0=a+0.5*h;for j=0:n-1qh=qh+fhs(x0+j*h);endt(i+1)=0.5*(t(i)+h*qh);qh=0;if i>=4for k=k:is(k)=(4*t(k+1)-t(k))/3;endfor m=m:k-1c(m)=(16*s(m+1)-s(m))/15;endfor z=z:m-1r(z)=(64*c(z+1)-c(z))/63;ende=abs(r1-r(z-1));r1=r(z-1);endif(e<0.00001)break;endendq=r(z-1);endglobal b a k m z e qh r1;k=0;m=0;z=0;e=1;qh=0;r1=0;b=3;a=1;[q]=romberg();disp(['输出结果为：r=',num2str(q)]);实验四C程序#include"stdio.h"#include"math.h"//已知量double x=0,h=0.0005,y[3]={0},p[4]={0.025,0.045,0.085,0.1}; double fx(double x,double y[3],int t);void main(){inti,t;double k[4],m[3];for(i=0;i<4;i++){for(x=0;x<=p[i];x=x+h){for(t=0;t<3;t++)m[t]=y[t];for(t=0;t<3;t++){k[0]=fx(x,m,t);m[t]=y[t]+0.5*h*k[0];k[1]=fx(x+0.5*h,m,t);m[t]=y[t]+0.5*h*k[1];k[2]=fx(x+0.5*h,m,t);m[t]=y[t]+h*k[2];k[3]=fx(x+h,m,t);//m[t]=y[t];y[t]=y[t]+h*(k[0]+2*k[1]+2*k[2]+k[3])/6;}}x=0;for(t=0;t<3;t++)printf("y%d(%lf)=%lf\t",t+1,p[i],y[t]);printf("\n");for(t=0;t<3;t++)y[t]=0;}}double fx(double x,double y[3],int t){double f;if(t==0)f=0*x+1;else if(t==1)f=0*x+y[2];elsef=0*x+1000-1000*y[1]-100*y[2];return f;}Matlab程序function[f]=fx(x,y)%UNTITLED7 Summary of this function goes here % Detailed explanation goes heref(1)=0*x+1;f(2)=0*x+y(3);f(3)=0*x+1000-1000*y(2)-100*y(3);endy=zeros(1,3);h=0.0005;x=0;p=[0.025,0.045,0.085,0.1];k=zeros(4,3);for i=1:4for j=1:fix(p(i)/h)k(1,:)=fx(x,y);k(2,:)=fx(x+0.5*h,y+0.5*h*k(1,:));k(3,:)=fx(x+0.5*h,y+0.5*h*k(2,:));k(4,:)=fx(x+h,y+h*k(3,:));y=y+h*(k(1,:)+2*k(2,:)+2*k(3,:)+k(4,:))/6;x=x+h;enddisp(['x=',num2str(p(i)),' ','y=',num2str(y)]);y=zeros(1,3);end实验五C程序#include"stdio.h"#include"math.h"#define N 9//已知量doubleA[N][N]={12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.0678 13,-2.031743,2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0 .037585,-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1 .784317,0.718719,1.563849,1.836742,-3.741856,0.431637,9.789356,-0.103458,-1.103456,0.238 417,1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2. 213474,3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00 001};doubleb[9]={2.1874369,33.992318,-25.173417,-0.84671695,1.784317,-86.612343,1.1101230,4. 719345,-5.6784392};//函数声明void change(inta,int b);//主函数void main(){ inti,j,q1,n,m,t;double tmp1,tmp2=0,x[9],tmp3;for(i=0;i<9;i++){tmp1=A[i][i];q1=i;for(j=i;j<9;j++)//{if(fabs(tmp1)<fabs(A[j][i])){ tmp1=A[j][i];q1=j;}}if(q1!=i)change(i,q1);for(n=i+1;n<9;n++){tmp3=A[i][i]/A[n][i];for(m=i+1;m<9;m++)A[n][m]=A[n][m]*tmp3-A[i][m];//需检查,为何主元用公式置零不行。
