Simpson 1/3 en Matlab
%regla de simpson 1/3
clear all; close all; clc
fun=input('Ingresa la funcin f(x) entre comillas: ');
f=inline(fun);
n=1;
while mod(n,2)~=0
n=input('Ingrese el nmero de subintervalos: ');
if mod(n,2)~=0
disp('El nmero de subintervalos debe ser par, pulse una tecla para continuar')
pause
end
end
a=input('Ingrese el lmite inferior de la integral: ');
b=input('Ingrese el lmite superior de la integral: ');
h=(b-a)/n;
sumai=0;
sumap=0;
for i=1:2:n-1
sumai=sumai+feval(f,h*i+a);
end
for i=2:2:n-2
sumap=sumap+feval(f,h*i+a);
end
int=(h/3)*(feval(f,a)+4*sumai+2*sumap+feval(f,b));
disp(['El resultado de la integral es ' num2str(int)])
%****************************************************************
%** Calculo de la integral por el
**
%** Metodo de Simpson de 1/3 y 3/8 para
**
%** multiples segmentos pares o impares
UdeG **
%**
Maestria en Electronica **
%** Ing. Jesus Norato Valencia
**
%** Materia: Metodos Numericos
**
%** Maestro: M.C. J.Gilberto Mateos
26/Nov/99 **
%****************************************************************
clear;
clc;
fprintf('Calculo de la integral por el metodo de Simpson de 1/3 y 3/8\n\n');
fprintf('Para cualquier cantidad de segmentos par o impar\n\n');
f=input('introduce la funcion:','s');
a=input('lime inferior:');
b=input('limite superior:');
c=input('numero de segmentos a dividir :');
h=(b-a)/c;
%*****************************************************************
%** En la siguiente seccion se aplica la regla de simpson 3/8
**
%** para los ultimos 3 segmentos (si es que la cantidad de
**
%** segmentos es impar de otra manera no entra a la sig. secc.) **
%*****************************************************************
qqqqq=0;
if (-1)^c==-1
b=b-(3*h);
c=c-3;
e=b+h;
ff=b+2*h;
g=b+3*h;
ee=g-b;
x=b;
q=eval(f);
x=e;
qq=eval(f);
x=ff;
qqq=eval(f);
x=g;
qqqq=eval(f);
qqqqq=(ee)*((q+3*qq+3*qqq+qqqq)/8);
end
%*****************************************************************
%** En la siguiente seccion se realiza la sumatoria
**
%** de cada una de las evaluaciones impares de la funcion
**
%*****************************************************************
z=0;
x=a;
for i=1:c;
if (-1)^i==1
k=eval(f);
z=z+k;
end
x=h*i;
end
%*****************************************************************
%** En la siguiente seccion se realiza la sumatoria
**
%** de cada una de las evaluaciones pares de la funcion
**
%*****************************************************************
zz=0;
x=a;
for i=2:c;
if (-1)^i==-1
k=eval(f);
zz=zz+k;
end
x=h*i;
end
%*****************************************************************
%** En la siguiente seccion se realiza la evaluacion de
**
%** el primero y el ultimo valor a evaluar en la funcion
**
%*****************************************************************
x=a;
if x==a
end
d=eval(f);
x=b;
if x==b
eee=eval(f);
end
%******************************************************************
%** una vez que se tienen los datos de areas bajo la curva
**
%** se realizan las operaciones directamente en la formula de
**
%** integracion por el metodo de simpson de 1/3 y se suma con
**
%** el valor de optenido de las areas de los ultimos 3 segmentos **
%******************************************************************
z=z*4;
v=zz*2;
z=z+v+d+eee;
z=z/(3*c);
z=z*(b-a);
z=z+qqqqq
fprintf('Resultado ');
Mtodo de Runge-Kutta
Programa en MatLab de Runge-Kutta de orden dos.
function f
fprintf('\n \tRESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO RUNGE-KUTTA DE ORDEN 4\n')
f=input('\n Ingrese la ecuacion diferencial dy/dx=\n','s');
x0=input('\n Ingrese el primer punto x0:\n');
x1=input('\n Ingrese el segundo punto x1:\n');
y0=input('\n Ingrese la condicion inicial y(x0):\n');
n=input('\n Ingrese el numero de pasos n:\n');
h=(x1-x0)/n;
xs=x0:h:x1;
fprintf('\n''it x0 y(x1)');
for i=1:n
it=i-1;
x0=xs(i);
x=x0;
y=y0;
k1=h*eval(f);
x=xs(i+1);
y=y0+k1;
k2=h*eval(f);
y0=y0+(k1+k2)/2;
fprintf('\n%2.0f%10.6f%10.6f\n',it,x0,y0);
end
fprintf('\n El punto aproximado y(x1) es = %8.6f\n',y0);
Programa de Runge-Kutta de orden cuatro.
function f
fprintf('\n \tRESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO RUNGE-KUTTA DE ORDEN 4\n')
f=input('\n Ingrese la ecuacion diferencial\n','s');
x0=input('\n Ingrese el primer punto x0:\n');
x1=input('\n Ingrese el segundo punto x1:\n');
y0=input('\n Ingrese la condicion inicial y(x0):\n');
n=input('\n Ingrese el numero de pasos n:\n');
h=(x1-x0)/n;
xs=x0:h:x1;
fprintf('\n''it x0 y(x1)');
for i=1:n
it=i-1;
x0=xs(i);
x=x0;
y=y0;
k1=h*eval(f);
x=x0+h/2;
y=y0+k1/2;
k2=h*eval(f);
x=x0+h/2;
y=y0+k2/2;
k3=h*eval(f);
x=x0+h;
y=y0+k3;
k4=h*eval(f);
y0=y0+(k1+2*k2+2*k3+k4)/6;
fprintf('\n%2.0f%10.6f%10.6f\n',it,x0,y0);
end
fprintf('\n El punto aproximado y(x1) es = %8.6f\n',y0);
Solucion
2. Respuesta
>> runge2
RESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO RUNGE-KUTTA DE ORDEN 4
Ingrese la ecuacion diferencial dy/dx=
4*exp(0.8*x)-0.5*y
Ingrese el primer punto x0:
0
Ingrese el segundo punto x1:
4
Ingrese la condicion inicial y(x0):
2
Ingrese el numero de pasos n:
4
'it x0 y(x1)
0 0.000000 6.701082
1 1.000000 16.319782
2 2.000000 37.199249
3 3.000000 83.337767
El punto aproximado y(x1) es = 83.337767
respuesta
>> runge4
RESOLUCION DE ECUACIONES DIFERENCIALES POR MEDIO RUNGE-KUTTA DE ORDEN 4
Ingrese la ecuacion diferencial
-2*x^3+12*x^2-20*x+8.5
Ingrese el primer punto x0:
0
Ingrese el segundo punto x1:
0.5
Ingrese la condicion inicial y(x0):
1
Ingrese el nmero de pasos n:
5
'it x0 y(x1)
0 0.000000 1.753950
1 0.100000 2.331200
2 0.200000 2.753950
3 0.300000 3.043200
4 0.400000 3.218750
function [T,Z] = rks4(Fn,a,b,Za,m)
%--------------------------------------------------------------------------%RKS4 Runge-Kutta solution for the system
% of D.E.'s Z' = F(t,Z) with Z(a) = Za.
% Sample call
% [T,Z] = rks4('f',a,b,ya,m)
% Inputs
% f name of the function
% a left endpoint of a<=t<=b
% b right endpoint of a<=t<=b
% Za vector of initial values
% m number of steps
% Return
% T parameter vector
% Z matrix of the vector solutions
%
%--------------------------------------------------------------------------h = (b - a)/m;
T = zeros(1,m+1);
Z = zeros(m+1,length(Za));
T(1) = a;
Z(1,:) = Za;
for j=1:m,
tj = T(j);
Zj = Z(j,:);
K1 = h*feval(Fn,tj,Zj);
K2 = h*feval(Fn,tj+h/2,Zj+K1/2);
K3 = h*feval(Fn,tj+h/2,Zj+K2/2);
K4 = h*feval(Fn,tj+h,Zj+K3);
Z(j+1,:) = Zj + (K1 + 2*K2 + 2*K3 + K4)/6;
T(j+1) = a + h*j;
end
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
Usage: [y t] = rk4(f,a,b,ya,n) or y = rk4(f,a,b,ya,n)
Runge-Kutta method of order 4 for initial value problems
Input:
f - Matlab inline function f(t,y)
a,b - interval
ya - initial condition
n - number of subintervals (panels)
Output:
y - computed solution
t - time steps
Examples:
[y t]=rk4(@myfunc,0,1,1,10);
here 'myfunc' is a user-defined function in M-file
y=rk4(inline('sin(y*t)','t','y'),0,1,1,10);
f=inline('sin(y(1))-cos(y(2))','t','y');
y=rk4(f,0,1,1,10);
function [y t] = rk4(f,a,b,ya,n)
h = (b - a) / n;
halfh = h / 2;
y(1,:) = ya;
t(1) = a;
h6 = h/6;
for i = 1 : n
t(i+1) = t(i) + h;
th2 = t(i) + halfh;
s1 = f(t(i), y(i,:));
s2 = f(th2, y(i,:) + halfh * s1);
s3 = f(th2, y(i,:) + halfh * s2);
s4 = f(t(i+1), y(i,:) + h * s3);
y(i+1,:) = y(i,:) + (s1 + s2+s2 + s3+s3 + s4) * h6;
end;
