Matlab6
Friday, 16 April 2021 1:56 PM Created a file
% demo How to open file to write
% File opening in a mode specified by permission. permission can be:
% r Open file for reading (default).
% w Open file, or create new file, for writing; discard existing contents, if any.
% a Open file, or create new file, for writing; append data to the end of the file.
% r+ open (do not create) file for reading and writing.
% w+ Open file, or create a new file, for reading and writing; discard existing contents, if any.
% a+ Open file, or create new file, for reading and writing; append data to the end of the file.
% Open a new file named IIEST1.txt at given location with write permission
fileID = fopen('/Users/prithadas/Documents/MATLAB/IIEST1.txt', 'w');
fprintf(fileID,'BESU is former name of IIEST, Shibpur, Howrah');
% Also check at folder location for the file created
fclose(fileID);
%type IIEST1.txt % Type the newly created file
File operation
% fileID = fopen(filename) opens the file, filename, for binary read access,
and returns
% an integer file identifier equal to or greater than 3.
% file identifiers is 0, 1, and 2 for standard input, standard output (the
screen), and standard error, respectively.
% If fopen cannot open the file, then fileID is -1
% Open the file, IIEST.txt, and obtain the file identifier.
% This txt file sld be where the the m file is.
fileID = fopen('IIEST1.txt');
%Pass the fileID to the fgetl function to read one line from the file.
tline = fgetl(fileID)
fclose(fileID); %close the file.
% if file is loacted in other place, absolute address can be used
fileID = fopen('/Users/prithadas/Documents/MATLAB/IIEST1.txt');
tline = fgetl(fileID)
fclose(fileID); %close the file.
IIEST1 text file
BESU is former BE College at Shibpur, Howrah
NOTE
File access type, specified as a character vector or a string scalar. You can open a file in binary mode or in
text mode. On UNIX systems, both translation modes have the same effect. To open a file in binary mode,
specify one of the following.
'r' Open file for reading.
'w' Open or create new file for writing. Discard existing contents, if any.
'a' Open or create new file for writing. Append data to the end of the file.
'r+' Open file for reading and writing.
'w+' Open or create new file for reading and writing. Discard existing contents, if any.
'a+' Open or create new file for reading and writing. Append data to the end of the file.
'A' Open file for appending without automatic flushing of the current output buffer.
'W' Open file for writing without automatic flushing of the current output buffer.
demo Lorenz eqn solved using ODE45
function X = lorenz(rho, sigma, beta, initV, T, eps)
% LORENZ Function generates the lorenz attractor of the prescribed values
% of parameters rho, sigma, beta
% Evaluates the right hand side of the Lorenz system
% x' = sigma*(y-x)
% y' = x*(rho - z) - y
% z' = x*y - beta*z
% parameters
sigma=10; rho=28; beta=2.6;
xst1 = sqrt(beta*(rho -1));
yst1= sqrt(beta*(rho -1));
zst1 = rho-1
xst2 =- sqrt(beta*(rho -1));
yst2= -sqrt(beta*(rho -1));
zst2 = rho-1
T = [0 25]; % Time span
initV = [0 1 1.05]% Intial values
[T,X] = ode45(@(T,X) F(T, X, sigma, rho, beta), T, initV);
figure('color','w');
plot3(X(:,1),X(:,2),X(:,3));
title('Lorenz Attractor','FontSize',10,'FontWeight','bold');
xlabel('X','FontSize',10,'FontWeight','bold');
xlabel('Y','FontSize',10,'FontWeight','bold');
xlabel('Z','FontSize',10,'FontWeight','bold');
set(gca,'LineWidth',2,'FontSize',10,'FontWeight','bold');
hold on;
scatter3(xst1,yst1,zst1, 50, 'r', 'filled');
scatter3(xst2,yst2,zst2, 50, 'r', 'filled');
%return
%end
function dx = F(T, X, sigma, rho, beta)
dx = zeros(3,1);
dx(1) = sigma*(X(2) - X(1));
dx(2) = X(1)*(rho - X(3)) - X(2);
function demoddepd()
% y'1(t)=y1(t-1) dx(3) = X(1)*X(2) - beta*X(3);
return
% y'2(t)=y1(t-1)+ y2(t-0.2)
% y'3(t)=y2(t)
end
% The history function for t0 is constant, y1(t)=y2(t)=y3(t)=1
end
% The time delays in the equations are only present in y terms, and the delays are % constants, so the equations form a system of constant delay
equations.
% to code the equations, delays, and history before calling the delay differential % equation (DDE) solver dde23
% First, create a vector to define the delays in the system of equations. This
% system has two different delays:
% A delay of 1 in the first component y1(t-1)
% A delay of 0.2 in the second component y2(t-0.2)
% dde23 accepts a vector argument for the delays, where each element is the
% constant delay for one component.
lags = [1 0.2];
% define the interval of integration [t0??tf]
tspan = [0 5];
% Solve the DDE using the dde23 solver
sol = dde23(@ddefun, lags, @history, tspan);
% Plot Solution
% The solution structure sol has the fields sol.x and sol.y that contain the
% internal time steps taken by the solver and corresponding solutions at those times.
% Plot the three solution components against time.
plot(sol.x,sol.y)
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2','y_3');
% Code Solution History
% Next, create a function to define the solution history. The solution history
% is the solution for times t<=t0
function s = history(t)
s = ones(3,1);
end
% Code Equation
% Now, create a function to code the equations. This function should have the
% signature dydt = ddefun(t,y,Z), where t is time (independent variable).
% y is the solution (dependent variable) Z(:,j) approximates the
% delay y(t- tau_j), where the constant delay tau_j is given by lags(j).
% These inputs are automatically passed to the function by the solver, but the
% variable names determine how you code the equations. In this case:
% Z(:,1)-> y1(t-1)
% Z(:,2)-> y2(t-0.2)
function dydt = ddefun(t,y,Z)
ylag1 = Z(:,1);
ylag2 = Z(:,2);
dydtpdesolve()
function = [ylag1(1);
% demo toylag1(1)+ylag2(2);
solve Partial Delay Differential Eqn (PDE)
y(2)];
% read documentation https://in.mathworks.com/help/matlab/math/solve-single-pde.html
endgiven as pi^2 del_u/del_t = del_^2u/del_x^2
% PDE
end
% The equation is defined on the interval 0<=x<=1 for times t>=0. At t=0,
% the solution satisfies the initial condition u(x,0)=sin(?x).
% V Also, at x=0 and x=1, the solution satisfies the boundary conditions
% u(0,t)=0,
% pi*e^-t + del_u/del_x(1,t)=0
% To solve this equation in MATLAB, you need to code the equation, the initial
% conditions, and the boundary conditions, then select a suitable solution mesh
% before calling the solver pdepe.
% Rewrite it in a form that the pdepe solver expects. The standard form that pdepe expects is
% c(x,t,u,del_u/del_x)del_u/del_t =x^-m del/del_x(x^m f(x,t,u,del_u/del_x))+s(x,t,u,?u/?x)
% Written in this form, the given PDE becomes
% pi^2 del_u/del_t = x^0 del_/del_x(x^0 del_u/del_x)+0
% With the equation in the proper form you can read off the relevant terms:
% m=0
% c(x,t,u,del_u/del_x)= pi^2
% f(x,t,u,del_u/del_x) = del_u/del_x
% s(x,t,u,del_u/del_x)= 0
% Now you can create a function to code the equation. The function should
% have the signature [c,f,s] = pdex1pde(x,t,u,dudx)
% For this problem, use a mesh with 20 equally spaced points in the spatial
% interval [0,1] and five values of t from the time interval [0,2].
% Select Solution Mesh
% Before solving the equation you need to specify the mesh points (t,x) at which
% you want pdepe to evaluate the solution. Specify the points as vectors t and x.
% The vectors t and x play different roles in the solver. In particular, the cost
% and accuracy of the solution depend strongly on the length of the vector x.
% However, the computation is much less sensitive to the values in the vector t.
x = linspace(0,1,20);
t = linspace(0,2,5);
% x is the independent spatial variable.
% t is the independent time variable.
% u is the dependent variable being differentiated with respect to x and t.
% dudx is the partial spatial derivative del_u/del_x.
% The outputs c, f, and s correspond to coefficients in the standard PDE equation form
% expected by pdepe. These coefficients are coded in terms of the input
% variables x, t, u, and dudx
% As a result, the equation in this example can be represented by the function:
% Finally, solve calling main PDE solver pdepe the equation using the symmetry m, the PDE equation, the initial
% conditions, the boundary conditions, and the meshes for x and t.
m = 0;
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% pdepe returns the solution in a 3-D array sol, where sol(i,j,k) approximates the
% kth component of the solution uk evaluated at t(i) and x(j). The size of sol is
% length(t)-by-length(x)-by-length(u0), since u0 specifies an initial condition for
% each solution component. For this problem, u has only one component, so
% sol is a 5-by-20 matrix, but in general you can extract the kth solution
% component with the command u = sol(:,:,k).
% Extract the first solution component from sol.
u = sol(:,:,1);
% Plot Solution: Create a surface plot of the solution.
surf(x,t,u)
title('Numerical solution computed with 20 mesh points')
xlabel('Distance x')
ylabel('Time t')
function [c,f,s] = pdex1pde(x,t,u,dudx)
c = pi^2;
f = dudx;
s = 0;
end
% Code Initial Condition
% Next, write a function that returns the initial condition. The initial condition
% is applied at the first time value tspan(1).
% The function should have the signature u0 = pdex1ic(x).
% The corresponding function here is
function u0 = pdex1ic(x)
u0 = sin(pi*x);
end
% Code Boundary Conditions
% Now, write a function that evaluates the boundary conditions. For problems
% posed on the interval a<= x <=b, the boundary conditions apply for all t and either
% x=a or x=b. The standard form for the boundary conditions expected by the solver is
% p(x,t,u)+q(x,t)f(x,t,u,del_u/del_x)=0
% Rewrite the boundary conditions in this standard form and read off the coefficient values.For x=0, the equation is
% u(0,t)=0 -> u+0.del_u/del_x =0
% The coefficients are:
% p(0,t,u)=u
% q(0,t)=0
% For x=1, the equation is
% pi*e^-t + del_u/del_x(1,t)= 0 -> pi*e^-t +1.del_u/del_x(1,t)= 0
% The coefficients are:
% p(1,t,u)=?e^-t
% q(1,t)=1
% Since the boundary condition function is expressed in terms of f(x,t,u,del_u/del_x), and
% this term is already defined in the main PDE function, you do not need to specify
% this piece of the equation in the boundary condition function. You need only specify
% the values of p(x,t,u) and q(x,t) at each boundary.
% The boundary function should use the function signature
%[pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t):
% The inputs xl and ul correspond to u and x for the left boundary.
% The inputs xr and ur correspond to u and x for the right boundary.
% t is the independent time variable.
% The outputs pl and ql correspond to p(x,t,u) and q(x,t)
% for the left boundary (x=0 for this problem).
% The outputs pr and qr correspond to p(x,t,u) and q(x,t) for the right boundary (x=1 for this problem).
% The boundary conditions in this example are represented by the function:
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;
end
end
