Matlab - Asistente tecnico

   
Vista:

Asistente tecnico

Publicado por Marcos (1 intervención) el 08/04/2009 15:20:20
Hola, estoy teniendo un problemita con la convergencia del Metodo de Newthon Rhapson para sistemas de ecuaciones lineales. Lo que trato de obtener es apartir de una ecuacion dendo no conosco 3 parametros, pero tengo muchos puntos como para comprobarlo. Los parametros a averiguar son k, a y b. D, E son constantes y dp, P y hp son variables de las cuales tengo varios puntos.
dp/D=k*(P*E*D^-2)^a * (hp/D)^b

A continuacion le copio el codigo:

clear,clc,cla,close all
format long
tol=0.01;
imax=200;
x=[0.1,0.1,0.1];
[xr,kp]=newtonsi(x,tol,imax)

-----------------------------------------------

function [xr,kp]=newtonsi(x,tol,imax)
format long
kp=1;
epi=1;
x1=x;
while norm(epi)>tol
x=x1;
fxn=f(x)
axn=jac(x)
epi=axnfxn';
x1=x-epi';
kp=kp+1
if kp>imax
disp('no converge')
break
end
end
xr=x1;

----------------------------------------------------

function y=f(x)
% funcion para utilizar con newtonsi.m
format long
E=200000;
p=429.005;
D=1.5875;
hp=0.04256;
dp=0.53642;

y(1)=x(1)*((p*E*D^(-2))^x(2))*(hp/D)^x(3)-dp/D

E=200000;
p=1030;
D=1.5875;
hp=0.0877;
dp=0.7596;

y(2)=x(1)*((p*E*D^(-2))^x(2))*(hp/D)^x(3)-dp/D

E=200000;
p=1900;
D=1.5875;
hp=0.18;
dp=1.10988;

y(3)=x(1)*((p*E*D^(-2))^x(2))*(hp/D)^x(3)-dp/D

----------------------------------------------------------

function df=jac(x)
% matriz jacobiana para usar con newtonsi.m
format long

E=200000;
p=429.005;
D=1.5875;
hp=0.04256;
dp=0.53642;

df(1,1)=(p*E/D^2)^x(2)*(hp/D)^x(3);
df(1,2)=x(1)*(p*E/D^2)^x(2)*log(p*E/D^2)*(hp/D)^x(3);
df(1,3)=x(1)*(p*E/D^2)^x(2)*(hp/D)^x(3)*log(hp/D);

E=200000;
p=1030;
D=1.5875;
hp=0.0877;
dp=0.7596;

df(2,1)=(p*E/D^2)^x(2)*(hp/D)^x(3);
df(2,2)=x(1)*(p*E/D^2)^x(2)*log(p*E/D^2)*(hp/D)^x(3);
df(2,3)=x(1)*(p*E/D^2)^x(2)*(hp/D)^x(3)*log(hp/D);

E=200000;
p=1900;
D=1.5875;
hp=0.18;
dp=1.10988;

df(3,1)=(p*E/D^2)^x(2)*(hp/D)^x(3);
df(3,2)=x(1)*(p*E/D^2)^x(2)*log(p*E/D^2)*(hp/D)^x(3);
df(3,3)=x(1)*(p*E/D^2)^x(2)*(hp/D)^x(3)*log(hp/D);
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder
información
Otras secciones de LWP con contenido de Matlab
- Código fuente de Matlab
- Cursos de Matlab
- Temas de Matlab
- Chat de Matlab
información
Códigos de Matlab
- Fuerzas de Coulomb
- TRIÁNGULO
- Elige una carta