Autor Tema: Misc

0 Usuarios y 1 Visitante están viendo este tema.

14 Enero, 2021, 01:05 am
Leído 396 veces

Ricardo Boza

  • $$\Large \color{#c88359}\pi\,\pi\,\pi\,\pi$$
  • Mensajes: 757
  • País: es
  • Karma: +1/-0
  • Sexo: Masculino

Código: (Matlab) [Seleccionar]
function fnomov(columnas)
% columnas=1,2,3 etc.
lon=8;
figure(1)
close(1)
figure(1)
hold on
axis equal
axis([-2 2 0 lon -4 1])
view(35,20)
n=20;
vec=linspace(-1,1,n);
for j=columnas:-1:1
    for i=1:n
        t=linspace(0,lon,50);
        fval1=evalua1(t,[vec(i),vec(j)]);
        plot3(vec(i)*ones(1,length(t)),t,vec(j)*ones(1,length(t))+fval1,'b')
        P=[vec(i),0,vec(j)];
        plot3(P(1),P(2),P(3),'bo','markerfacecolor','b')
        plot3([P(1) P(1)],[0,lon],[P(3) P(3)],'k')
        pause(0.1)
    end
end
pause(0.9)
v=linspace(35,450,150);
for k=1:150
    view(v(k),20)
    pause(0.1)
end
pause(0.3)
v=linspace(0,20,50);
for k=1:50
    view(90,v(51-k))
    pause(0.1)
end
pause(1)
view(35,20)

    function[y1]=evalua1(t0,y0)
        [~,y]=ode45(@fun,t0,y0);
        y1=y(:,1);
    end

    function[dydt]=fun(~,y)
        y1=-y(1)+3*y(2);
        y2=-y(2);
        dydt=[y1; y2];
    end
end

[attachment id=0 msg=460163]

Código: (Matlab) [Seleccionar]
function fcno
y0=linspace(0,2,20);
t0=linspace(0,6,601);
figure(1)
close(1)
figure(1)
options=odeset('RelTol',1e-5);
subplot(1,3,1)
hold on
for k=1:20
    [t,y]=ode45(@fun,t0,y0(k),options);
    plot(t,y,'b')
end
title('ode')
subplot(1,3,2)
hold on
t=t0;
for k=1:20
    f2=y0(k)*exp(t)./(1+(exp(t)-1)*y0(k));
    plot(t,f2,'r')
end
title('solution')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lon=601;
t=linspace(0,6,lon);
h=6/(lon-1);
y0=linspace(0,2,20);
subplot(1,3,3)
hold on
for m=1:20
    y=zeros(1,lon);
    y(1)=y0(m);
    for k=1:lon-1
        y(k+1)=y(k)+h*y(k)*(1-y(k));
    end
    plot(t,y,'g')
end
plot([0 6],[0 0],'m--','linewidth',2)
plot([0 6],[1 1],'m--','linewidth',2)
title('my ode')

    function [yp]=fun(~,y)
        yp=y*(1-y);
    end
end

[attachment id=1 msg=460163]

Código: (Matlab) [Seleccionar]
a=1.03;
h=(a^200-a^20)/9;
ypoints=a^(20):h:a^200;
xpoints=log(ypoints)./log(a);

figure(1)
close(1)
figure(1)
hold on
axis([0 220 -1 1])
plot([20 200],[0 0],'b-','linewidth',3)
z=0.05;
for k=1:10
plot([xpoints(k) xpoints(k)],z*[-1 1],'b','linewidth',2)
end
text(20-1,-1.5*z,'20')
text(200-2,-1.5*z,'200')


Código: (Matlab) [Seleccionar]
function fourier1(n)
% Copy and paste this on Command Window:
% for n=1:15, fourier1(n), end
% Or this:
% fourier1(30)
figure(1)
clf
T=2;
lon=200;
a_0=2/T*integral(@(x)x,-T/2,T/2);
a=zeros(1,n);
b=zeros(1,n);
for k=1:n
    a(k)=2/T*integral(@(x)x.*cos(2*pi*k/T*x),-T/2,T/2);
    b(k)=2/T*integral(@(x)x.*sin(2*pi*k/T*x),-T/2,T/2);
end

x=linspace(-T/2,T/2,lon);
y=a_0/2*ones(1,lon);

for k=1:n
    y=y+b(k)*sin(2*pi*k/T*x);
end

figure(1)
hold on
% plot([-1 1],[-1 1],'k')
subplot(1,3,1)
hold on
axis equal
axis([-1 1 -1 1])
plot([-1 1],[-1 1],'k')
plot(x,y,'b')
title(['n=' num2str(n)])

for k=1:n
    c(k)=(-1)^(k+1)*T/(pi*k);
end

z=zeros(1,lon);
for k=1:n
    z=z+c(k)*sin(2*pi*k/T*x);
end
subplot(1,3,2)
hold on
axis equal
axis([-1 1 -1 1])
plot([-1 1],[-1 1],'k')
plot(x,z,'r')
title(['n=' num2str(n)])

subplot(1,3,3)
hold on
axis equal
axis([-1 1 -1 1])
plot([-1 1],[-1 1],'k')
plot(x,y,'b')
plot(x,z,'k.')
title(['n=' num2str(n)])

pause(1)
end



Código: (Matlab) [Seleccionar]
figure(1)
close(1)
figure(1)
T=2;
lon=200;
x=linspace(-T/2,T/2,lon);
for n=1:10
    clf
    a_0=2/T*integral(@(x)x.^2,-T/2,T/2);
    a=zeros(1,n);
    b=zeros(1,n);
    for k=1:n
        a(k)=2/T*integral(@(x)x.^2.*cos(2*pi*k/T*x),-T/2,T/2);
        b(k)=2/T*integral(@(x)x.^2.*sin(2*pi*k/T*x),-T/2,T/2);
    end
   
    y=a_0/2*ones(1,lon);
    for k=1:n
        y=y+a(k)*cos(2*pi*k/T*x)+b(k)*sin(2*pi*k/T*x);
    end
   
    figure(1)
    hold on
    axis([-1 1 -0.2 1])
    plot([-1 1],[0 0],'k--')
    plot([0 0],[-0.2 1],'k--')
    plot(x,x.^2,'r')
    plot(x,y,'b')
    title(['n=' num2str(n)])
    pause(1)
end


Código: (Matlab) [Seleccionar]
function fourier_num(order)
% Copy and paste in Command Window:
% fourier_num(50)

T=2;
lon=500;
x=linspace(-T/2,T/2,lon);
figure(1)
close(1)
figure(1)
for n=order:order
    a_0=2/T*integral(@(x)fun1(x),-T/2,T/2,'ArrayValued',true);
    a=zeros(1,n);
    b=zeros(1,n);
    for k=1:n
        a(k)=2/T*integral(@(x)fun1cos(x,k,T),-T/2,T/2,'ArrayValued',true);
        b(k)=2/T*integral(@(x)fun1sin(x,k,T),-T/2,T/2,'ArrayValued',true);
    end
   
    y=a_0/2*ones(1,lon);
    for k=1:n
        y=y+a(k)*cos(2*pi*k/T*x)+b(k)*sin(2*pi*k/T*x);
    end
   
    clf
    hold on
    axis equal
    axis([-1.5 1.5 -1.5 1.5])
    plot([-1.5 1.5],[0 0],'k--')
    plot([0 0],[-1.5 1.5],'k--')
    % plot(x,x.^2,'r')
    plot(x,y,'b')
    title(['n=' num2str(n)])
    pause(0.2)
end

    function[y]=fun1(x)
        if x<0
            y=-1;
        else
            y=1;
        end
    end

    function[y]=fun1cos(x,k,T)
        if x<0
            y=-cos(2*pi*k/T*x);
        else
            y=cos(2*pi*k/T*x);
        end
    end

    function[y]=fun1sin(x,k,T)
        if x<0
            y=-sin(2*pi*k/T*x);
        else
            y=sin(2*pi*k/T*x);
        end
    end

end



Código: (Matlab) [Seleccionar]
function fourier_num_pause(order)
% Copy and paste in Command Window:
% fourier_num_pause(30)

T=2;
lon=500;
x=linspace(-T/2,T/2,lon);
figure(1)
close(1)
figure(1)
for n=1:order
    a_0=2/T*integral(@(x)fun1(x),-T/2,T/2,'ArrayValued',true);
    a=zeros(1,n);
    b=zeros(1,n);
    for k=1:n
        a(k)=2/T*integral(@(x)fun1cos(x,k,T),-T/2,T/2,'ArrayValued',true);
        b(k)=2/T*integral(@(x)fun1sin(x,k,T),-T/2,T/2,'ArrayValued',true);
    end
   
    y=a_0/2*ones(1,lon);
    for k=1:n
        y=y+a(k)*cos(2*pi*k/T*x)+b(k)*sin(2*pi*k/T*x);
    end
   
    clf
    hold on
    axis equal
    axis([-1.5 1.5 -1.5 1.5])
    plot([-1.5 1.5],[0 0],'k--')
    plot([0 0],[-1.5 1.5],'k--')
    % plot(x,x.^2,'r')
    plot(x,y,'b')
    title(['n=' num2str(n)])
    pause(0.2)
end

    function[y]=fun1(x)
        if x<0
            y=-1;
        else
            y=1;
        end
    end

    function[y]=fun1cos(x,k,T)
        if x<0
            y=-cos(2*pi*k/T*x);
        else
            y=cos(2*pi*k/T*x);
        end
    end

    function[y]=fun1sin(x,k,T)
        if x<0
            y=-sin(2*pi*k/T*x);
        else
            y=sin(2*pi*k/T*x);
        end
    end

end



Código: (Matlab) [Seleccionar]
clc
lon=500;
T=2;
figure(1)
close(1)
figure(1)
%%%%%%%%%%%%%
% Define f
%%%%%%%%%%%%%
for n=1:15
    syms f(x)
    % f(x)=x;
    % f(x)=piecewise(x<0, -1, 0<x, 1);
    % f(x)=x^2;
    f(x)=x^3;
    syms fcos(x,k)
    fcos(x,k)=f(x)*cos(2*pi*k/T*x);
    syms fsin(x,k)
    fsin(x,k)=f(x)*sin(2*pi*k/T*x);
    a_0=2/T*int(f(x),-T/2,T/2);
    a=sym(zeros(1,n));
    b=sym(zeros(1,n));
    for k=1:n
        a(k)=2/T*int(fcos(x,k),-T/2,T/2);
        b(k)=2/T*int(fsin(x,k),-T/2,T/2);
    end
   
    y=a_0/2;
    for k=1:n
        y=y+a(k)*cos(2*pi*k/T*x)+b(k)*sin(2*pi*k/T*x);
    end
   
    x=linspace(-T/2,T/2,100);
    clf
    hold on
    axis equal
    axis([-1.5 1.5 -1.5 1.5])
    plot([-1.5 1.5],[0 0],'k--')
    plot([0 0],[-1.5 1.5],'k--')
    plot(x,subs(f),'k')
    plot(x,subs(y),'b')
    title(['n=' num2str(n)])
    pause(0.1)
end
disp('y=')
pretty(y)



Código: (Matlab) [Seleccionar]
clc
lon=500;
T=2;
figure(1)
close(1)
figure(1)
%%%%%%%%%%%%%
% Define f
%%%%%%%%%%%%%
for m=1:8 % <------------------ Modify this parameter
    for n=1:15
        syms f(x)
        if m==1
            f(x)=x;
        elseif m==2
            f(x)=piecewise(x<0, -1, 0<x, 1);
        elseif m==3
            f(x)=x^2;
        elseif m==4
            f(x)=x^3;
        elseif m==5
            f(x)=piecewise(x<0,-x, 0<x, x);
        elseif m==6
            f(x)=piecewise(x<0,cos(pi*x),0<x,-cos(pi*x));
        elseif m==7
            f(x)=piecewise(x<0,cos(pi*x),0<x, cos(2*pi*x));
        elseif m==8
            f(x)=piecewise(x<0,sin(pi*x),0<x<0.5,-2*x+1,0.5<x,2*x-1);
        end
        syms fcos(x,k)
        fcos(x,k)=f(x)*cos(2*pi*k/T*x);
        syms fsin(x,k)
        fsin(x,k)=f(x)*sin(2*pi*k/T*x);
        a_0=2/T*int(f(x),-T/2,T/2);
        a=sym(zeros(1,n));
        b=sym(zeros(1,n));
        for k=1:n
            a(k)=2/T*int(fcos(x,k),-T/2,T/2);
            b(k)=2/T*int(fsin(x,k),-T/2,T/2);
        end
       
        y=a_0/2;
        for k=1:n
            y=y+a(k)*cos(2*pi*k/T*x)+b(k)*sin(2*pi*k/T*x);
        end
       
        x=linspace(-T/2,T/2,100);
        clf
        hold on
        axis equal
        axis([-1.5 1.5 -1.5 1.5])
        plot([-1.5 1.5],[0 0],'k--')
        plot([0 0],[-1.5 1.5],'k--')
        plot(x,subs(f),'k')
        plot(x,subs(y),'b')
        title(['n=' num2str(n)])
        pause(0.1)
    end
end
text(1,-1,'end')



 
[attachment id=2 msg=460163]

Código: (Matlab) [Seleccionar]
function exer1
n=20;
tmax=6;
t=linspace(0,tmax,20);
y0=linspace(-1,1,n);
cte=1-1./y0.^2;
figure(1)
close(1)
figure(1)
for k=1:n
    [t,z]=ode45(@(t,y)fun(t,y),t,y0(k));
    subplot(1,3,1)
    hold on
    axis equal
    axis([0 tmax -1 1])
    plot(t,z,'bo-')
    y=1./sqrt((1-cte(k)*exp(2*t)));
    subplot(1,3,2)
    hold on
    axis equal
    axis([0 tmax -1 1])
    plot(t,y,'ko-')
    plot(t,-y,'ko-')
    subplot(1,3,3)
    hold on
    axis equal
    axis([0 tmax -1 1])
    plot(t,z,'b')
    plot(t,y,'k')
    plot(t,-y,'k')
end
figure(1)
subplot(1,3,1)
title('ode45')
subplot(1,3,2)
title('solution')
subplot(1,3,3)
title('comparison')
    function [dydt]=fun(~,y)
        dydt=y*(y^2-1);
    end
end