正文
sor函数python python中的solve函数
提示：扫一扫查出行【扫一扫了解最新限行尾号】
复制提示
请c/c++高手看一下关于sor程序的问题
错误不在sor，在这里：
for(j=0;jm;j++)
{
for(l=0;ln;l++)
{
a[j][l]=rand();
b[j][l]=rand();
c[j][l]=rand();
d[j][l]=rand();
e[j][l]=rand();
f[j][l]=rand();
u[j][l]=0.0;
}
}
怀疑你的dmatrix函数写的有问题，代码写得很乱，不知道在干什么。估计是内存分配有问题。我这里的运行时错误是访问了非法地址0。环境是VS 2005.
我把你的程序稍微修改了一下，可以看得比较清楚错误的地方。
#include stdlib.h
#include stdio.h
//#include "nrutil.h"
#include time.h
#include math.h
#include stddef.h
#define MAXITS 1000
#define NR_END 1
#define FREE_ARG char*
#define EPS 1.0e-5
void nrerror(char error_text[]);
void sor(double **a,double **b,double **c,double **d,double **e,double **f,double **u,int jmax,double rjac)
{
int ipass,j,jsw,l,lsw,n;
double anorm,anormf=0.0,omega=1.0,resid;
printf("enter sor\n");
for(j=2;jjmax;j++)
for(l=2;ljmax;l++)
anormf+=fabs(f[j][l]);
for(n=1;n=MAXITS;n++)
{
anorm=0.0;
jsw=1;
for(ipass=1;ipass=2;ipass++)
{
lsw=jsw;
for(j=2;jjmax;j++)
{
for(l=lsw+1;ljmax;l+=2)
{
resid=a[j][l]*u[j+1][l]+b[j][l]*u[j-1][l]+c[j][l]*u[j][l+1]+d[j][l]*u[j][l-1]+e[j][l]*u[j][l]-f[j][l];
anorm+=fabs(resid);
u[j][l]-=omega*resid/e[j][l];
}
lsw=3-lsw;
}
jsw=3-jsw;
omega=(n==1ipass==1?1.0/(1.0-0.5*rjac*rjac):1.0/(1.0-0.25*rjac*rjac*omega));
}
if(anormEPS*anormf)
return;
}
nrerror("MAXITS exceeded");
}
void nrerror(char error_text[])
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
double **dmatrix(long nrl, long nrh, long ncl, long nch)
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i=nrh;i++) m[i]=m[i-1]+ncol;
printf("dmatrix succeeded!\n");
return m;
}
main()
{
double **a,**b,**c,**d,**e,**f,**u,rjac;
int j,l,jmax;
double pi=3.1415926;
long m=10,n=10;
a=dmatrix(1,m,1,n);
b=dmatrix(1,m,1,n);
c=dmatrix(1,m,1,n);
d=dmatrix(1,m,1,n);
e=dmatrix(1,m,1,n);
f=dmatrix(1,m,1,n);
u=dmatrix(1,m,1,n);
printf("%x, %x, %x, %x, %x, %x, %x\n", (int)a, (int)b, (int)c, (int)d, (int)e, (int)f, (int)u);
scanf("%d",jmax);
rjac=cos(pi/jmax);
srand((unsigned)time(NULL));
printf("init begin...\n");
for(j=0;jm;j++)
{
for(l=0;ln;l++)
{
a[j][l]=rand();
b[j][l]=rand();
c[j][l]=rand();
d[j][l]=rand();
e[j][l]=rand();
f[j][l]=rand();
u[j][l]=0.0;
}
}
printf("init finished!\n");
sor(a,b,c,d,e,f,u,jmax,rjac);
for(j=0;j10;j++)
{
for(l=0;l10;l++)
{
printf("%f\t",u[j][l]);
}
printf("\n");
}
getch();
}
相关问答
Q1: python类的题目求帮忙看下
class Person(object):
name = "张三"
def __init__(self):
self.__age = 18
def set_age(self, age):
self.__age = age
#使用装饰器@classmethod修饰getName()方法
@classmethod
def getName(cls):
return cls.name
#使用类名访问name属性
print(Person().name)
#生成对象p1
p1 = Person()
#使用p1访问name属性
print(p1.name)
#给p1动态绑定属性country
p1.country = 'China'
#获取name
print(p1.getName())
Q2: python hashlib.sha512 为什么加密结果和shadow不一样
pythonsor函数python的话直接使用crypt就行计算sor函数python的时候要加上salt值。
cat /etc/shadow | greb root
得到密文。 密文由3部分组成，以”$”分隔，第一部分为ID，第二部分为盐值，第三部分为加密密文
密文示例：
$6$D0xsORq3b7GGsbYv$7L8myKFcJf1FqnvtzD.Zcbx9PeHnwGH1YJm/xDeEKmABPIxzUC7mQmEiIEHoHypXeJ969Tfp4c47WQgjkT0YH0
$id$salt$encrypted
采用方法 盐值 加密密文
$6$D0xsORq3b7GGsbYv #这一段就是盐值，在加密的时候连同前面的ID一起加密
import crypt
print crypt.crpyt("password","$6$D0xsORq3b7GGsbYv")
就可以得到shadow中的密文了。
其实shadow中的密文并不是直接由sha512得出的，还有一些其他操作,例如使用base64加密了salt和hash，可以看下这里man page
Q3: SOR迭代法 求MATLAB程序
迭代法 matlab实现代码如下
function [x,n] = jacobi(A,b,x0,eps,varargin)
if nargin ==3
eps = 1.0e-6;
M = 200;
elseif nargin3
disp('输入参数数目不足3个');
return
elseif nargin ==5
M = varargin{1};
end
D = diag(diag(A)); %%求Asor函数python的对角矩阵
L = -tril(A,-1); %%求Asor函数python的下三角矩阵
U = -triu(A,1); %%求A的上三角矩阵
B = D\(L+U);
f = D\b;
x = B*x0+f;
n = 1;%迭代次数
while norm(x-x0)=eps
x0 = x;
x = B*x0+f
n = n+1;
if(n=M)
disp('Warning:迭代次数太多,可能不收敛!')
return;
end
end
运行效果如下sor函数python：
扩展资料：
迭代法的收敛性判别
收敛性判别条件
SOR迭代法收敛的充分必要条件是ρ(λω)1sor函数python，ρ(λω)与松弛因子ω有关。ρ(λω)与ω的关系以及SOR方法收敛的条件有如下定理。
定理1：(Kahan)对任意的A
sor函数python，设其对角元皆非零，则对所有实数ω，有：ρ(λω)≥ ω-1。
推论：如果解Ax=b的SOR方法收敛，则有ω-11，即0ω2。
定理2：(Ostrowski-Reich)设A
，A对称正定，且0ω2，则解Ax=b的SOR方法收敛。
参考资料来源：百度百科-逐次超松驰迭代法
Q4: 如何用MATLAB编写的拉格朗日插值算法的程序、二阶龙格-库塔方法的程序和SOR迭代法的程序
拉格朗日function y=lagrange(x0,y0,x)
n=length(x0);m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end SOR迭代法的Matlab程序
function [x]=SOR_iterative(A,b)
% 用SOR迭代求解线性方程组,矩阵A是方阵
x0=zeros(1,length(b)); % 赋初值
tol=10^(-2); % 给定误差界
N=1000; % 给定最大迭代次数
[n,n]=size(A); % 确定矩阵A的阶
w=1; % 给定松弛因子
k=1;
% 迭代过程
while k=N
x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1);
for i=2:n
x(i)=(1-w)*x0(i)+w*(b(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')/A(i,i);
end
if max(abs(x-x0))=tol
fid = fopen('SOR_iter_result.txt', 'wt');
fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n');
fprintf(fid,'迭代次数: %d次\n\n',k);
fprintf(fid,'x的值\n\n');
fprintf(fid, '%12.8f \n', x);
break;
end
k=k+1;
x0=x;
end
if k==N+1
fid = fopen('SOR_iter_result.txt', 'wt');
fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n');
fprintf(fid,'迭代次数: %d次\n\n',k);
fprintf(fid,'超过最大迭代次数，求解失败！');
fclose(fid);
end Matlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。
(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'）
Yn+1=Yn+h*f(Xn,Yn)
另外根据微分中值定理，存在0t1,使得
Yn+1=Yn+h*f(Xn+th,Y(Xn+th))
这里K＝f(Xn+th,Y(Xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。
利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h^5)的四阶龙格库塔公式：
K1＝f(Xn,Yn);
K2=f(Xn+h/2,Yn+(h/2)*K1);
K3=f(Xn+h/2,Yn+(h/2)*K2);
K4=f(Xn+h,Yn+h*K3);
Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6);
所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。
仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数：
function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）
n=floor((b-a)/h);%求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值，可以是向量，但是要注意维数
for ii=1:n
x(ii+1)=x(ii)+h;
k1=ufunc(x(ii),y(:,ii));
k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);
k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);
k4=ufunc(x(ii)+h,y(:,ii)+h*k3);
y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
%按照龙格库塔方法进行数值求解
end
调用的子函数以及其调用语句：
function dy=test_fun(x,y)
dy = zeros(3,1);%初始化列向量
dy(1) = y(2) * y(3);
dy(2) = -y(1) + y(3);
dy(3) = -0.51 * y(1) * y(2);
对该微分方程组用ode45和自编的龙格库塔函数进行比较，调用如下：
[T,F] = ode45(@test_fun,[0 15],[1 1 3]);
subplot(121)
plot(T,F)%Matlab自带的ode45函数效果
title('ode45函数效果')
[T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%测试时改变test_fun的函数维数，别忘记改变初始值的维数
subplot(122)
plot(T1,F1)%自编的龙格库塔函数效果
title('自编的 龙格库塔函数')
关于sor函数python和python中的solve函数的介绍到此就结束了，不知道你从中找到你需要的信息了吗 ？如果你还想了解更多这方面的信息，记得收藏关注本站。
