### Autor Tema: Misc

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

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

#### Ricardo Boza

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

Código: (Matlab) [Seleccionar]
function fnomov(columnas)% columnas=1,2,3 etc.lon=8;figure(1)close(1)figure(1)hold onaxis equalaxis([-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)    endendpause(0.9)v=linspace(35,450,150);for k=1:150    view(v(k),20)    pause(0.1)endpause(0.3)v=linspace(0,20,50);for k=1:50    view(90,v(51-k))    pause(0.1)endpause(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];    endend
[attachment id=0 msg=460163]

Código: (Matlab) [Seleccionar]
function fcnoy0=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 onfor k=1:20    [t,y]=ode45(@fun,t0,y0(k),options);    plot(t,y,'b')endtitle('ode')subplot(1,3,2)hold ont=t0;for k=1:20    f2=y0(k)*exp(t)./(1+(exp(t)-1)*y0(k));    plot(t,f2,'r')endtitle('solution')%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%lon=601;t=linspace(0,6,lon);h=6/(lon-1);y0=linspace(0,2,20);subplot(1,3,3)hold onfor 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')endplot([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);    endend
[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 onaxis([0 220 -1 1])plot([20 200],[0 0],'b-','linewidth',3)z=0.05;for k=1:10plot([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)clfT=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);endx=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);endfigure(1)hold on% plot([-1 1],[-1 1],'k')subplot(1,3,1)hold onaxis equalaxis([-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);endz=zeros(1,lon);for k=1:n    z=z+c(k)*sin(2*pi*k/T*x);endsubplot(1,3,2)hold onaxis equalaxis([-1 1 -1 1])plot([-1 1],[-1 1],'k')plot(x,z,'r')title(['n=' num2str(n)])subplot(1,3,3)hold onaxis equalaxis([-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    endend
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    endend
Código: (Matlab) [Seleccionar]
clclon=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)enddisp('y=')pretty(y)
Código: (Matlab) [Seleccionar]
clclon=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)    endendtext(1,-1,'end')

[attachment id=2 msg=460163]

Código: (Matlab) [Seleccionar]
function exer1n=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')endfigure(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);    endend