数值分析五个题目的C语言及Matlab程序
数值分析MATLAB编程-推荐下载
题目三：设������(������) = ������������，在[-1,1]上用 Legendre 多项式作������(������)的 3 次最佳平方逼近多项式。 >> a=-1;b=1; >> n=3; >> syms x; >> fx=exp(x); >> exp(x); >> for k=3 F=x^k*fx; d(k+1)=int(F,a,b); end >> d=reshape(d,n+1,1); >> H=hilb(n+1); >Fra bibliotek a=H\d a=
对全部高中资料试卷电气设备，在安装过程中以及安装结束后进行高中资料试卷调整试验；通电检查所有设备高中资料电试力卷保相护互装作置用调与试相技互术通关，1系电过，力管根保线据护敷生高设产中技工资术艺料0不高试仅中卷可资配以料置解试技决卷术吊要是顶求指层，机配对组置电在不气进规设行范备继高进电中行保资空护料载高试与中卷带资问负料题荷试2下卷2，高总而中体且资配可料置保试时障卷，各调需类控要管试在路验最习；大题对限到设度位备内。进来在行确管调保路整机敷使组设其高过在中程正资1常料中工试，况卷要下安加与全强过，看度并22工且22作尽22下可22都能22可地护以缩1关正小于常故管工障路作高高；中中对资资于料料继试试电卷卷保破连护坏接进范管行围口整，处核或理对者高定对中值某资，些料审异试核常卷与高弯校中扁对资度图料固纸试定，卷盒编工位写况置复进.杂行保设自护备动层与处防装理腐置，跨高尤接中其地资要线料避弯试免曲卷错半调误径试高标方中高案资等，料，编试要5写、卷求重电保技要气护术设设装交备备置底4高调、动。中试电作管资高气，线料中课并敷3试资件且、设卷料中拒管技试试调绝路术验卷试动敷中方技作设包案术，技含以来术线及避槽系免、统不管启必架动要等方高多案中项；资方对料式整试，套卷为启突解动然决过停高程机中中。语高因文中此电资，气料电课试力件卷高中电中管气资壁设料薄备试、进卷接行保口调护不试装严工置等作调问并试题且技，进术合行，理过要利关求用运电管行力线高保敷中护设资装技料置术试做。卷到线技准缆术确敷指灵设导活原。。则对对：于于在调差分试动线过保盒程护处中装，高置当中高不资中同料资电试料压卷试回技卷路术调交问试叉题技时，术，作是应为指采调发用试电金人机属员一隔，变板需压进要器行在组隔事在开前发处掌生理握内；图部同纸故一资障线料时槽、，内设需，备要强制进电造行回厂外路家部须出电同具源时高高切中中断资资习料料题试试电卷卷源试切，验除线报从缆告而敷与采设相用完关高毕技中，术资要资料进料试行，卷检并主查且要和了保检解护测现装处场置理设。备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况，然后根据规范与规程规定，制定设备调试高中资料试卷方案。
数值分析中求解线性方程组的MATLAB程序(6种)
数值分析中求解线性方程组的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程序
数值分析（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一、准备工作我们需要明确2021 数学建模C 题的具体要求和背景。
在这次比赛中，题目要求参赛选手利用给定的数据集和背景信息，通过建立数学模型，并编写相关代码进行模拟和求解。
针对不同的场景和要求，我们需要结合 MATLAB 软件的特点和优势，从数据处理、模型建立到结果分析，全面展现我们的解题思路和方法。
二、数据处理及可视化分析在开始编写 MATLAB 代码之前，我们需要对给定的数据集进行分析和处理。
具体来说，可以通过 MATLAB 的数据导入工具，将原始数据导入到 MATLAB 评台中，并进行数据预处理、清洗和可视化分析。
在处理数据的过程中，可以选用 MATLAB 中丰富的数据处理函数和绘图工具，展现数据的分布特点、趋势规律以及异常情况。
三、模型建立和求解在数据处理的基础上，我们将重点转向模型的建立和求解。
在C 题中，可能涉及到不同的数学模型，比如优化模型、拟合模型、动态模型等。
我们可以逐步引入相关的 MATLAB 函数和工具箱，比如优化工具箱、符号计算工具箱等，来构建模型并进行求解。
在代码的编写过程中，需要考虑到代码的优化和可读性，以便审阅人员和评审团队能够清晰地理解我们的解题思路和方法。
四、结果分析与讨论我们需要对模型的求解结果进行分析和讨论。
通过 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。
数值分析 算法C语言程序
一、拉格朗日插值#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班姓名：甄彦福。
数值分析的MATLAB程序
列主元法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')。
数值分析五个题目的C语言及MATLAB程序
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() {
break; end end 实验二 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,
} 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 here for i=1:9
h(i)=x(i+1)-x(i); a(i)=(fx(i)-fx(i+1))/(x(i)-x(i+1)); j=i-1;
end end function [q] =ytdhs(m,h,p,fx,x) %UNTITLED14 Summary of this function goes here % Detailed explanation goes here for i=1:10
(完整word版)同济大学数值分析matlab编程题汇编
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 ≤≤上求解初值问题。
数值分析matlab代码
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方法求下列微分方程的解。
数值分析大作业(牛顿下山法,拉格朗日法,切比雪夫法)及Matlab程序
课程设计课程名称：数值分析设计题目：学号：姓名：完成时间：2014.11.18题目一： 解线性方程组的直接法 设方程组Ax b =，其中250002511125555111x x x x x x A x x x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦， 矩阵中10.1(0,1,,5)k x k k =+=，b 由相应的矩阵元素计算，使解向量(1,1,,1)T x =。
(1) A 不变，对b 的元素6b 加一个扰动410-，求解方程组；(2) b 不变，对A 的元素22a 和66a 分别加一个扰动610-，求解方程组； (3) 对上述两种扰动方程组的解做误差分析。
一.数学原理：本计算采用直接法中的列主元高斯消元法，高斯列主元消元法原理如下： 1、设有n 元线性方程组如下：1111n n nn a a a a ⎛⎫ ⎪ ⎪ ⎪⎝⎭1nx x ⎛⎫ ⎪ ⎪ ⎪⎝⎭＝1nb b ⎛⎫ ⎪ ⎪ ⎪⎝⎭2、第一步：如果a11!=0, 令l i1= ai1/a11, I= 2,3,……,n用（-li1）乘第一个方程加到第i 个方程上,得同解方程组：a (1)11 a (1)12 . . . a (1)1nx 1 b (1)1 a (1)21 a (1)22 . . . a (1)2n x 2 b (1)2 . . . . . . . = . a (1)n-11 a (1)n-12 . . a (1)n-1n x n-1 b (1)n-1 a (1)n1 a (1)n2 . . . a (1)nn x n b (1)n简记为：A (2) x = b (2) 其中a (2)ij = a (1)ij – l i1 * a (1)1j , I ,j = 2,3,..,nb (2)I = b (1)I – l i1 * b (1)1 , I = 2,3,...,n 第二步：如果a (2)22 != 0,令l i2= a （2）i2/a （2）22, I= 3,……,n依据同样的原理，对矩阵进行化间（省略），依次下去，直到完成！最后，得到上三角方程组：a(1)11 a(1)12. . . a(1)1nx1b(1)10 a(1)22 . . . a(1)2nx2b(1)2. . . . . . . = .0 0 . . a(n-1)n-1n xn-1b(n-1)n-10 0 . . . a(n)nn xnb(n)n简记为：A(n) x = b(n)最后从方程组的最后一个方程进行回代求解为：Xn = b(n) / a(n)nnXi = ( b(k)k- ∑ a(k)kj x j ) / a(k)kk二．解题过程：1.由题中所给条件可求出b。
数值分析与MATLAB课程 实验程序
第二次实验内容：首先下载10.22.74.73网站下数值计算课程“实验内容”栏目里的第二章进行学习和实践（其中“二、矩阵的特殊参数运算；三、矩阵的分解；四、矩阵中的一些特殊处理函数”部分暂时不看），然后完成：1） 定义-4pi<=x<=4pi ，y1=sinx/x, y2=2sin(x+0.4pi)/(x+0.4pi)，在同一坐标系中画出上面两条曲线，两条曲线用不同的颜色，不同的连线标记，不同的顶点标记来加以区分。
（提示：用plot 指令，可用help plot 学习MATLAB 的在线帮助）2） 输入如下两个矩阵 A 和 B ，对矩阵 A 和 B 作关系运算，标识出两矩阵中元素相等的位置，元素值不等的位置，并标识出矩阵 A 中所有小于 0 的元素。
3) 对上面的矩阵 A 和 B 作逻辑“ 或”、“ 与”运算，并标识出矩阵 B 中所有大于 2并小于 5 的元素位置。
4) 利用矩阵的寻址方式，取出上面矩阵A 中第二行元素作为新矩阵的第一列，B 矩阵的第3列元素作为为新矩阵的第二列元素，取出矩阵A 的第一行，第三行元素作为新矩阵的第三列和第四列元素，1. x=-4*pi:0.1:4*pi; y1=sin(x)./x;y2=2*sin(x+0.4*pi)./(x+0.4*pi);//数组所以是用点除。
应该是数组和数组的运算用点 plot(x,y1,'-black')//作函数y1黑色- hold on//hold 住plot(x,y2,'red+') //作函数y2红色+ title('曲线')//标注 xlabel('x') ylabel('y')hold off//不用hold 了 2. A=[1 2 3;-2 1 3;-3 2 1]; B=[1 4 3;3 2 8;5 2 3];A==B//这是逻辑运算，符合1，不符合0 A~=B A<03. A=[1 2 3;-2 1 3;-3 2 1]; B=[1 4 3;3 2 8;5 2 3]; A|B A&B(B>2)&(B<5)4.A=[1 2 3;-2 1 3;-3 2 1]; B=[1 4 3;3 2 8;5 2 3];C=[(A(2,:))' B(:,3) (A(1,:))' (A(3,:))'] 运行结果：123213321A ⎡⎤⎢⎥=-⎢⎥⎢⎥-⎣⎦143328523B ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦1.2. ans =1 0 1 0 0 0 0 1 0 ans =0 1 0 1 1 1 1 0 1 ans =0 0 0 1 0 0 1 0 0 3. ans =1 1 1 1 1 1 1 1 1 ans =1 1 1 1 1 1 1 1 1 ans =0 1 1 1 0 0 0 0 1 4. C =-2 3 1 -3 1 8 2 2 3 3 3 1 第二次实验内容： 对于非线性方程2()sin 102x x π---=（1） 编写M –File 函数用二分法求出其在区间［0，0.5］和［2.5，3］内的根，要求函数的最大循环次次为1000次，两个精度均为10-4。
matlab简单程序运用-数值分析篇
1：三次样条插值法h0=figure('toolbar','none',...'position',[200 50 350 450],...'name','实例86');h1=axes('parent',h0,...'position',[0.10 0.45 0.8 0.5],...'visible','off');x=0:0.2:2*pi;y=sin(x);plot(x,y)b1=uicontrol('parent',h0,...'units','points',...'tag','b1',...'style','pushbutton',...'string','三次样条插值',...'backgroundcolor',[0.75 0.75 0.75],...'position',[20 60 70 20],...'callback',[...'y=0,',...'sy=0,',...'strn1=get(e2,''string'');,',...'n1=str2num(strn1);,',...'strn2=get(e3,''string'');,',...'n2=str2num(strn2);,',...'x=n1:0.2:n2;,',...'i=1;,',...'for t=n1:0.2:n2,',...'y(i)=sin(t);,',...'sy(i)=san(t,n1,n2);,',...'i=i+1;,',...'end,',...'plot(x,y,''b*'',x,sy,''r-''),',...'axis([0 7 -1.5 1.5]),',...'legend(''sin(x)'',''N-Hermite插值'')']); b2=uicontrol('parent',h0,...'units','points',...'tag','b2',...'style','pushbutton',...'string','误差比较',...'backgroundcolor',[0.75 0.75 0.75],...'position',[170 60 70 20],...'callback',[...'strdn1=get(e2,''string'');,',...'n1=str2num(strdn1);,',...'strdn2=get(e3,''string'');,',...'n2=str2num(strdn2);,',...'strdn=get(e1,''string'');,',...'dn=str2num(strdn);,',...'dd=abs(sin(dn)-san(dn,n1,n2));,',...'msgbox([''误差为：'',num2str(dd)],''计算结果'')']); e1=uicontrol('parent',h0,...'units','points',...'tag','e1',...'style','edit',...'fontsize',12,...'string','1.20',...'horizontalalignment','right',...'backgroundcolor',[1 1 1],...'position',[200 100 40 20]);t1=uicontrol('parent',h0,...'units','points',...'tag','t1',...'style','text',...'string','误差点：',...'fontsize',12,...'backgroundcolor',[0.75 0.75 0.75],...'position',[160 100 40 20]);e2=uicontrol('parent',h0,...'units','points',...'tag','e2',...'style','edit',...'fontsize',12,...'string','1.00',...'horizontalalignment','right',...'backgroundcolor',[1 1 1],...'position',[20 85 40 20]);t2=uicontrol('parent',h0,...'units','points',...'tag','t2',...'style','text',...'string','第一节点：',...'fontsize',12,...'backgroundcolor',[0.75 0.75 0.75],...'position',[15 105 50 20]);e3=uicontrol('parent',h0,...'units','points',...'tag','e3',...'fontsize',12,...'string','3.00',...'horizontalalignment','right',...'backgroundcolor',[1 1 1],...'position',[100 85 40 20]);t3=uicontrol('parent',h0,...'units','points',...'tag','t3',...'style','text',...'string','第二节点：',...'fontsize',12,...'backgroundcolor',[0.75 0.75 0.75],...'position',[95 105 50 20]);b3=uicontrol('parent',h0,...'units','points',...'tag','b3',...'style','pushbutton',...'string','关闭',...'backgroundcolor',[0.75 0.75 0.75],...'position',[100 20 60 20],...'callback','close');2：NEWTON插值h0=figure('toolbar','none',...'position',[200 50 350 450],...'name','实例87');h1=axes('parent',h0,...'position',[0.10 0.45 0.8 0.5],...'visible','off');x=0:0.2:2*pi;y=sin(x);plot(x,y)b1=uicontrol('parent',h0,...'units','points',...'tag','b1',...'style','pushbutton',...'string','牛顿插值',...'backgroundcolor',[0.75 0.75 0.75],...'position',[20 60 70 20],...'callback',[...'strn=get(e1,''string'');,',...'n=str2num(strn);,',...'x=0:0.2:2*pi;,',...'for t=0:0.2:2*pi,',...'y(i)=sin(t);,',...'ynt(i)=newton(t,n);,',...'i=i+1;,',...'end,',...'plot(x,y,''b*'',x,ynt,''r-''),',...'axis([0 7 -1.5 1.5]),',...'legend(''sin(x)'',''牛顿插值'')']);b2=uicontrol('parent',h0,...'units','points',...'tag','b2',...'style','pushbutton',...'string','误差比较',...'backgroundcolor',[0.75 0.75 0.75],...'position',[170 60 70 20],...'callback',[...'strn=get(e1,''string'');,',...'n=str2num(strn);,',...'strdn=get(e2,''string'');,',...'dn=str2num(strdn);,',...'dd=abs(sin(dn)-newton(dn,n));,',...'msgbox([''误差为：'',num2str(dd)],''计算结果'')']); e1=uicontrol('parent',h0,...'units','points',...'tag','e1',...'style','edit',...'fontsize',12,...'string','5',...'horizontalalignment','right',...'backgroundcolor',[1 1 1],...'position',[50 100 40 20]);e2=uicontrol('parent',h0,...'units','points',...'tag','e2',...'style','edit',...'fontsize',12,...'string','1.20',...'horizontalalignment','right',...'backgroundcolor',[1 1 1],...'position',[200 100 40 20]);t1=uicontrol('parent',h0,...'units','points',...'tag','t1',...'string','节点数：（<6）',...'fontsize',12,...'backgroundcolor',[0.75 0.75 0.75],...'position',[10 100 40 30]);t2=uicontrol('parent',h0,...'units','points',...'tag','t2',...'style','text',...'string','误差点：',...'fontsize',12,...'backgroundcolor',[0.75 0.75 0.75],...'position',[160 100 40 20]);b3=uicontrol('parent',h0,...'units','points',...'tag','b3',...'style','pushbutton',...'string','关闭',...'backgroundcolor',[0.75 0.75 0.75],...'position',[100 20 60 20],...'callback','close');3mewton形式的hermite插值h0=figure('toolbar','none',...'position',[200 50 350 450],...'name','实例89');h1=axes('parent',h0,...'position',[0.10 0.45 0.8 0.5],...'visible','off');x=0:0.2:2*pi;y=sin(x);plot(x,y)b1=uicontrol('parent',h0,...'units','points',...'tag','b1',...'style','pushbutton',...'string','N-Hermite插值',...'backgroundcolor',[0.75 0.75 0.75],...'position',[20 60 70 20],...'callback',[...'strn1=get(e2,''string'');,',...'n1=str2num(strn1);,',...'strn2=get(e3,''string'');,',...'n2=str2num(strn2);,',...'x=0:0.2:2*pi;,',...'i=1;,',...'for t=0:0.2:2*pi,',...'y(i)=sin(t);,',...'ynh(i)=nhermite(t,n1,n2);,',...'i=i+1;,',...'end,',...'plot(x,y,''b*'',x,ynh,''r-''),',...'axis([0 7 -1.5 1.5]),',...'legend(''sin(x)'',''N-Hermite插值'')']);b2=uicontrol('parent',h0,...'units','points',...'tag','b2',...'style','pushbutton',...'string','误差比较',...'backgroundcolor',[0.75 0.75 0.75],...'position',[170 60 70 20],...'callback',[...'strdn1=get(e2,''string'');,',...'n1=str2num(strdn1);,',...'strdn2=get(e3,''string'');,',...'n2=str2num(strdn2);,',...'strdn=get(e1,''string'');,',...'dn=str2num(strdn);,',...'dd=abs(sin(dn)-nhermite(dn,n1,n2));,',...'msgbox([''误差为：'',num2str(dd)],''计算结果'')']); e1=uicontrol('parent',h0,...'units','points',...'tag','e1',...'style','edit',...'fontsize',12,...'string','1.20',...'horizontalalignment','right',...'backgroundcolor',[1 1 1],...'position',[200 100 40 20]);t1=uicontrol('parent',h0,...'units','points',...'tag','t1',...'style','text',...'string','误差点：',...'fontsize',12,...'backgroundcolor',[0.75 0.75 0.75],...'position',[160 100 40 20]);e2=uicontrol('parent',h0,...'units','points',...'tag','e2',...'style','edit',...'fontsize',12,...'string','1.00',...'horizontalalignment','right',...'backgroundcolor',[1 1 1],...'position',[20 85 40 20]);t2=uicontrol('parent',h0,...'units','points',...'tag','t2',...'style','text',...'string','第一节点：',...'fontsize',12,...'backgroundcolor',[0.75 0.75 0.75],...'position',[15 105 50 20]);e3=uicontrol('parent',h0,...'units','points',...'tag','e3',...'style','edit',...'fontsize',12,...'string','3.00',...'horizontalalignment','right',...'backgroundcolor',[1 1 1],...'position',[100 85 40 20]);t3=uicontrol('parent',h0,...'units','points',...'tag','t3',...'style','text',...'string','第二节点：',...'fontsize',12,...'backgroundcolor',[0.75 0.75 0.75],...'position',[95 105 50 20]);b3=uicontrol('parent',h0,...'units','points',...'tag','b3',...'style','pushbutton',...'string','关闭',...'backgroundcolor',[0.75 0.75 0.75],...'position',[100 20 60 20],...'callback','close');。
《数值分析》上机实验报告
数值分析上机实验报告《数值分析》上机实验报告1.用Newton 法求方程 X 7-X 4+14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。
1.1 理论依据：设函数在有限区间[a ，b]上二阶导数存在，且满足条件{}αϕ上的惟一解在区间平方收敛于方程所生的迭代序列迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)(')()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20)()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f ab c f x f b a x f b f x f k k k k k k ==-==∈≤-≠>+令)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3225333647>⋅''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f故以1.9为起点⎪⎩⎪⎨⎧='-=+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。
当前后两个的差<=ε时，就认为求出了近似的根。
本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。
1.2 C语言程序原代码：#include<stdio.h>#include<math.h>main(){double x2,f,f1;double x1=1.9; //取初值为1.9do{x2=x1;f=pow(x2,7)-28*pow(x2,4)+14;f1=7*pow(x2,6)-4*28*pow(x2,3);x1=x2-f/f1;}while(fabs(x1-x2)>=0.00001||x1<0.1); //限制循环次数printf("计算结果：x=%f\n",x1);}1.3 运行结果：1.4 MATLAB上机程序function y=Newton(f,df,x0,eps,M)d=0;for k=1:Mif feval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))<=epsd=1;breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey= '奇异'endfunction y=df(x)y=7*x^6-28*4*x^3;Endfunction y=f(x)y=x^7-28*x^4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 问题讨论：1.使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围内得到比较理想的计算结果。
数值分析在生活中的应用举例及Matlab实现
一、最小二乘法，用MATLAB实现1. 数值实例下面给定的是乌鲁木齐最近1个月早晨7：00左右（新疆时间）的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。
下面用MATLAB编程对上述数据进行最小二乘拟合。
下面用MATLAB编程对上述数据进行最小二乘拟合2、程序代码x=[1:1:30];y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];a1=polyfit(x,y,3) %三次多项式拟合%a2= polyfit(x,y,9) %九次多项式拟合%a3= polyfit(x,y,15) %十五次多项式拟合%b1=polyval(a1,x)b2=polyval(a2,x)b3=polyval(a3,x)r1= sum((y-b1).^2) %三次多项式误差平方和%r2= sum((y-b2).^2) %九次次多项式误差平方和%r3= sum((y-b3).^2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1, 'r') %用红色线画出x,b1图像%hold onplot(x,b2, 'g') %用绿色线画出x,b2图像%hold onplot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%3、数值结果不同次数多项式拟合误差平方和为：r1=67.6659r2=20.1060r3=3.7952r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。
4、拟合曲线如下图二、 线性方程组的求解（ 高斯-塞德尔迭代算法 ）1、实例： 求解线性方程组（见书P233页）⎪⎪⎩⎪⎪⎨⎧=++=-+=+-3612363311420238321321321x x x x x x x x x 记A x=b, 其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=363320,,12361114238321b x A x x x任取初始值()()Tx0000=，进行迭代。
数值分析matlab数值试验
实验一：误差传播及算法稳定性实验1.21、试验程序：function charpt1_2% 误差传播及算法稳定性实验clc;clear all;promps={'请选择递推关系式，若选E1=1/e,En=1-nEn-1,请输入1，若选EN=0，En-1=(1-En)/n,请输入2：'};I=1;while Iresult=inputdlg(promps,'charpt1_2',1,{'1'});Nb=str2num(char(result));if ((Nb~=1)|(Nb~=2))I=0;endend%%%%%%%%%%%%%%%%%I=1;while Iresult=inputdlg('请输入递推步数 n>=1:','charpt1_2',1,{'10'});steps=str2num(char(result));if (steps>0)&(steps==fix(steps)) %% 如果steps大于0且为整数I=0;endend%%%%%%%%%%%%%%%%%result=inputdlg('请输入计算中所采用的有效数字位数n:','charpt1_2',1,{'5'});Sd=str2num(char(result));format long %% 设置显示精度result=zeros(1,steps); %% 存储计算结果err=result; %% 存储计算的绝对误差值func=result; %% 存储用quadl计算的近似值%%%%%%%%%%%%%%%%%%% 用quadl计算积分近似值for n=1:stepsfun=@(x) x.^n.*exp(x-1);func(n)= quadl(fun,0,1);end%%%%%%%%%%%%%%%%%%% 用自定义算法计算if(Nb==1)digits(Sd);result(1)=subs(vpa(1/exp(1)));for n=2:stepsresult(n)=subs(vpa(1-n*result(n-1)));enderr=abs(result-func);elseif(Nb==2)digits(Sd);result(steps)=0;for n=(steps-1):-1:1result(n)=subs(vpa((1-result(n+1))/(n+1)));enderr=abs(result-func);end%%%%%%%%%%%%%%%%%%% 输出结果数值及图像clf;disp('库函数计算值:');disp(sprintf('%e ',func));disp('递推值:');disp(sprintf('%e ',result));disp('误差值:');disp(sprintf('%e ',err));if(Nb==1)plot([1:steps],result,'-rs',[1:steps],func,':k*',[1:steps],err,'-.bo' );elseif(Nb==2)plot([steps:-1:1],result,'-rs',[steps:-1:1],func,':k*',[steps:-1:1],e rr,'-.bo');endxlabel('第n步');ylabel('计算值');legend('自定义算法结果','库函数计算结果','误差值');grid on2、试验结果：选择递推关系式1，递推步数为10，有效数字为5位，计算结果如下：库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-0011.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002递推值:3.678800e-001 2.642400e-001 2.072800e-001 1.708800e-001 1.456000e-0011.264000e-001 1.152000e-001 7.840000e-0022.944000e-001 -1.944000e+000误差值:5.588280e-007 1.117662e-006 3.352927e-006 1.341222e-0056.705713e-005 4.023702e-004 2.816427e-003 2.253226e-002 2.027877e-001 2.027877e+00012345678910第n 步计算值选择递推关系式2，递推步数为10，有效数字为5位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678800e-001 2.642400e-001 2.072800e-001 1.708900e-001 1.455300e-001 1.267900e-001 1.125000e-001 1.000000e-001 1.000000e-001 0.000000e+000 误差值:5.588280e-007 1.117662e-006 3.352927e-006 3.412224e-006 2.942873e-006 1.237016e-005 1.164270e-004 9.322618e-004 8.387707e-003 8.387707e-002第n 步计算值选择递推关系式1，递推步数为10，有效数字为6位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678790e-001 2.642420e-001 2.072740e-001 1.709040e-001 1.454800e-001 1.271200e-001 1.101600e-001 1.187200e-001 -6.848000e-002 1.684800e+000 误差值:4.411720e-007 8.823378e-007 2.647073e-006 1.058778e-0055.294287e-005 3.176298e-004 2.223573e-003 1.778774e-002 1.600923e-001 1.600923e+00012345678910第n 步计算值选择递推关系式2，递推步数为10，有效数字为6位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678800e-001 2.642410e-001 2.072770e-001 1.708930e-001 1.455360e-001 1.267860e-001 1.125000e-001 1.000000e-001 1.000000e-001 0.000000e+000 误差值:5.588280e-007 1.176622e-007 3.529274e-007 4.122239e-007 3.057127e-006 1.637016e-005 1.164270e-004 9.322618e-004 8.387707e-003 8.387707e-002第n 步计算值3、结果分析：很明显第二种递推式结果要比第一种好，式1在第七步后有明显误差，而式2在第三步后基本与近似解一致。
(完整word版)matlab数值分析例题
1、 在MATLAB 中用Jacobi 迭代法讨论线性方程组,1231231234748212515x 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 —10 0 —10 0 0B=inv（D）*（L+U）% B为雅可比迭代矩阵B =0 0.2500 —0。
25000.5000 0 0.12500。
4000 —0.2000 0〉〉r=eigs(B,1）%B的谱半径r =0。
3347 〈1Jacobi迭代法收敛。
(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。
00004.00003。
0000k =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）>=epsx0=x;x=B＊x0+f；k=k+1;if（k〉=max1)disp(’迭代超过300次，方程组可能不收敛’）；return；endend2、设有某实验数据如下：（1)在MATLAB中作图观察离散点的结构,用多项式拟合的方法拟合一个合适的多项式函数；（2）在MATLAB中作出离散点和拟合曲线图。
数值分析编程题c语言汇总
数值分析实习报告- 1 -数值分析实习报告上机实习题一一、题目：b与已知A12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.0317 432.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.0101 031.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585A=-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.7841370.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.2384171.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.21 34743.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.00001b={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345 ,-5.6784392}1．用Household变换，把A化为三对角阵B（并打印B）。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
本文档包含上一个文档中的五个数值分析实验题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];//需检查,为何主元用公式置零不行。