diff options
| -rw-r--r-- | Functions/Heun/Heun.m | 44 | ||||
| -rw-r--r-- | Functions/days/README.md | 11 | ||||
| -rw-r--r-- | Functions/days/days.m | 22 | ||||
| -rw-r--r-- | Functions/falsePosition/falsePosition.m | 45 | ||||
| -rw-r--r-- | Functions/luFactor/luFactor.m | 87 | ||||
| -rw-r--r-- | Functions/simpson/README.md | 0 | ||||
| -rw-r--r-- | Functions/simpson/simpson.m | 39 | ||||
| -rw-r--r-- | Functions/specialMatrix/specialMatrix.m | 51 | ||||
| -rw-r--r-- | README.md | 15 | ||||
| -rw-r--r-- | Scripts/hw11/hw11.m | 29 | ||||
| -rw-r--r-- | Scripts/hw2/hw2.m | 60 | ||||
| -rw-r--r-- | Scripts/hw3/hw3.m | 37 | ||||
| -rw-r--r-- | Test.txt | 1 |
13 files changed, 439 insertions, 2 deletions
diff --git a/Functions/Heun/Heun.m b/Functions/Heun/Heun.m new file mode 100644 index 0000000..5812dfe --- /dev/null +++ b/Functions/Heun/Heun.m @@ -0,0 +1,44 @@ +function [t,y] = Heun(dydt,tspan,y0,h,es,maxit) +% [t,y] = Heun(dydt,tspan,y0,h): uses the heun method to integrate an ODE +% inputs: +% dydt = the differential equation of interest (must be anonymous +% function) +% tspan = [ti,tf] the initial and final values of the independent +% variable +% y0 = the initial value of the dependent variable +% h = step size +% es = stopping criterion (%), optional (default = 0.001) +% maxit = maximum iterations of corrector, optional (default = 50) +% outputs: +% t = vector of independent variable values +% y = vector of solution for dependent variable + +%Error checking +if nargin<4, error('At least 4 aruments required.');end +if nargin<5|isempty(es), es = 0.001;end +if nargin<6|isempty(maxit), maxit = 50;end +if length(tspan)~=2, error('Argument tspan must have a length of 2');end + +%Trims tspan according to step size (h) +t=tspan(1):h:tspan(2); +%if rem(tspan(2),h)~=0 +if t(length(t))~=tspan(2) + t=[t,tspan(end)]; +end + +y=ones(1,length(t))*y0; +h=diff(t); + +for i=1:length(t)-1 + S=dydt(t(i),y(i)); + ynew=y(i)+S*h; + for j=1:maxit + yold=ynew; + ynew=y(i)+h(i)/2*(S+dydt(t(i+1),yold)); + ea = abs((ynew-yold)/ynew)*100; + if ea <= es, break, end + end + y(i+1)=ynew; +end +plot(t,y) +end
\ No newline at end of file diff --git a/Functions/days/README.md b/Functions/days/README.md new file mode 100644 index 0000000..8a11db2 --- /dev/null +++ b/Functions/days/README.md @@ -0,0 +1,11 @@ +#days.m +Function to count the total days elapsed in a year according to a given date. +The syntax of the function is `days(<months>, <days>, <leap>)` where `months` is a integer (1-12). `days` is an integer (1-31) and `leap` + +## Input +months - month number (1-12). Example: `8` represents August. +days - day number of the month. +leap - indicates if the year is a leapyear or a regular year. `0` for regular and `1` for leap year. + +## Output +`nd` - number of days elapsed in the year.
\ No newline at end of file diff --git a/Functions/days/days.m b/Functions/days/days.m new file mode 100644 index 0000000..87163b7 --- /dev/null +++ b/Functions/days/days.m @@ -0,0 +1,22 @@ +function nd = days(mo, da, leap) +%days - Counts days elapsed in the year. +% +% days(<months>, <days>, <leap>) determines the days elapsed in a year +% based on the current date. This will include the current day. +% Months (1-12), +% Example: August 28, no leap year +% days(8,28,0) + +daysPeM=[0 31 59 90 120 151 181 212 243 273 304 334]; % cummulative days prior each month. + +nd = daysPeM(mo) + da; + +if leap == 1 & mo >= 3 + nd = nd + 1; +end +end + + +% MECH 105: Homework 4 - Part 1 +% Author: Christian Kolset +% Date: 5.21.21
\ No newline at end of file diff --git a/Functions/falsePosition/falsePosition.m b/Functions/falsePosition/falsePosition.m new file mode 100644 index 0000000..d2cb776 --- /dev/null +++ b/Functions/falsePosition/falsePosition.m @@ -0,0 +1,45 @@ +function [root, fx, ea, iter] = falsePosition(func, xl, xu, es, maxit, varargin) +%falsePosition finds the root of a function using false position method +% Syntax: +% [root, fx, ea, iter] = falsePosition(func, xl, xu, es, maxit, varargin) +% +% Input: +% func = function being evaluated +% xl = lower guess +% xu = upper guess +% es = desired relative error (default 0.0001%) +% maxit = maximum number of iterations (default 200) +% varargin, ...= any additional parameters used by the function +% Output: +% root = estimated root location +% fx = function evaluated at root location +% ea = approximated relative error (%) +% iter = number of iterations performed + +if nargin < 3, error('Too few arguments: at least 3 needed'), end +if nargin < 4 | isempty(es), es = 0.0001; end +if nargin < 5 | isempty(maxit), maxit = 200; end + +iter=0; +ea=100; +xr= xl; + +while(1) + xrold= xr; + iter = iter + 1; + xr=xu-func(xu)*(xl-xu)/(func(xl)-func(xu)); + if xr~=0 + ea = abs((xr-xrold)/xr)*100; + end + if func(xl)*func(xr)<0 + xl=xr; + else + xl=xr; + end + if iter == maxit | ea <= es, break, end +end +root = xr; +fx=func(xr,varargin{:}); +ea=ea; +iter=iter; +end
\ No newline at end of file diff --git a/Functions/luFactor/luFactor.m b/Functions/luFactor/luFactor.m new file mode 100644 index 0000000..4c853dd --- /dev/null +++ b/Functions/luFactor/luFactor.m @@ -0,0 +1,87 @@ +function [L, U, P] = luFactor(A) +% luFactor(A) +% LU decomposition with pivoting +% inputs: +% A = coefficient matrix +% outputs: +% L = lower triangular matrix +% U = upper triangular matrix +% P = the permutation matrix + + +[m,n] = size(A); +if m~=n, error('Enter a square matrix.'); end +nb = n + 1; +L=zeros(m,m); +U=zeros(m,m); +P=eye(m); + +%Partial pivoting +for j = 1:n - 1 + for p = j+1:m + if (abs(A(j,j)) < abs(A(p,j))) + A([j p],:) = A([p j],:); + P([j p],:) = P([p j],:); + end + end +end +% [big,i] = max(abs(A(j:n,j))); +% ipr = i + j - 1; +% if ipr ~= j +% A([j,ipr],:) = A([ipr,j],:); +% end +%end + +for i=1:m + % Computing L + for k=1:i-1 + L(i,k)=A(i,k); + for j=1:k-1 + L(i,k)= L(i,k)-L(i,j)*U(j,k); + end + L(i,k) = L(i,k)/U(k,k); + end + + % Computing U + for k=i:m + U(i,k) = A(i,k); + for j=1:i-1 + U(i,k)= U(i,k)-L(i,j)*U(j,k); + end + end +end + +for i=1:m + L(i,i)=1; +end + +%% Checks +if P*A ~= L*U, error('LU factorization failed to be computed correctly.'); end + +% Checks output size +if size(L) ~= size(A) | size(U) ~= size(A) | size(P) ~= size(A) + error('Incorrect output size') +end + +% Checks if matrices are triangular +%for i=1:m +% for j=1:n +% error('Incorrect output: LU decomposition are not triangular') +%end + +% Checks the permutation matrix P is valid +%P_str = string(P) +%for k = m +% num1=count(P(:,1:n),1) +% if num1 ~=1 +% error('Incorrect P matrix output. Too many 1's on the same line') +% end +% num0=count(P(:,1:n),0) +% if num0 ~= n-1 +% error('Incorrect P matric output. Wrong number of 0's') +% end +%end + +%} + +end
\ No newline at end of file diff --git a/Functions/simpson/README.md b/Functions/simpson/README.md new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/Functions/simpson/README.md diff --git a/Functions/simpson/simpson.m b/Functions/simpson/simpson.m new file mode 100644 index 0000000..acd4a76 --- /dev/null +++ b/Functions/simpson/simpson.m @@ -0,0 +1,39 @@ +function [I] = Simpson(x, y) +% Numerical evaluation of integral by Simpson's 1/3 Rule +% Inputs +% x = the vector of equally spaced independent variable +% y = the vector of function values with respect to x +% Outputs: +% I = the numerical integral calculated + +%Error cheching +if length(x)~=length(y), error('The vectors need to be of equal lengths'),end +if length(x)==1, error('Vectors must be at least 2 data points long'),end +if range(x(2:end)-x(1:end-1))~=0, error('The x vector needs to be equally spaced'), end +%Error Errorchecking END + +lower_bound=min(x); +interval=length(x)-1; +if rem(interval,2)==0 + trap_rule=0; + upper_bound=max(x); +else + warning('Odd number of intervals. Applying Trapezoidal rule on last interval') + trap_rule=1; + upper_bound=x(end-1); +end + +if length(x)==2 + I=(upper_bound-lower_bound)*(y(1)+y(2))/2; +else + h=(upper_bound-lower_bound)/interval; + SumOdd=sum(y(1:2:end))-y(end-1); + SumEven=sum(y(2:2:end))-y(end-2); + I=(upper_bound-lower_bound)*(y(1)+4*SumOdd+2*SumEven+y(end))/(3*interval); +end + +if trap_rule==1 & length(x)~=2 + I_trap=(max(x)-upper_bound)*(y(upper_bound)+y(max(x)))/2; + I=I+I_trap; +end +end
\ No newline at end of file diff --git a/Functions/specialMatrix/specialMatrix.m b/Functions/specialMatrix/specialMatrix.m new file mode 100644 index 0000000..d6c31f8 --- /dev/null +++ b/Functions/specialMatrix/specialMatrix.m @@ -0,0 +1,51 @@ +function [A] = specialMatrix(n,m) +% This function should return a matrix A as described in the problem statement +% Inputs n is the number of rows, and m the number of columns +% It is recomended to first create the matrxix A of the correct size, filling it with zeros to start with is not a bad choice + +%-------------------------------------------- + +if nargin ~= 2 + error('Error: Please enter two arguments.') +end +if (n|m) <= 0 + error('Error: Index error. Arguments must be positive integers or logical values.') +end + +A = zeros(n,m); %Creates a "blank" n x m matrix full of zeros + +% Now the real challenge is to fill in the correct values of A + +A(1,:) = 1:m; %Lables the first row with column numbers +A(:,1) = 1:n; %Lables the first column with row numbers + +for j = 2:n + for k = 2:m + A(j,k)=A(j-1,k)+A(j,k-1); % fills each element of the matrix with the sum of the left and above element. + end +end +end + +%--------------------------------------------- + +%{ +if nargin(specialMatrix) ~= 2 + error('Error: Please enter two arguments.') +elseif nargin(specialMatrix) <= 0 + error('Error: Index errer. Argument must be positive integers or logical values.') +end + +if n<2 + A = [1:m]; %If matrix is 1 x m then make a matrix with 1 row. +else + A = [1:m;1:m]; %Create row 1 +end + +for j = 3:n %Loop to make + A(j,:) = sum(A); % +end % +A(:,1)=(1:n); + +end +% Things beyond here are outside of your functions +%}
\ No newline at end of file @@ -1,2 +1,15 @@ # mech105-algorithms -Algorithms I created in MECH-105 (Summer 2021) +Algorithms I created in MECH-105 (2021-Summer-Term). Consists of MATLAB scripts and functions + +##Function: +*specialMatrix.m +*falsePosition.m +*luFactor.m +*Heun.m +*simpson.m +*days.m + +##Scripts: +*hw2.m +*hw3.m +*hw11.m
\ No newline at end of file diff --git a/Scripts/hw11/hw11.m b/Scripts/hw11/hw11.m new file mode 100644 index 0000000..312e816 --- /dev/null +++ b/Scripts/hw11/hw11.m @@ -0,0 +1,29 @@ +% Define problem constants +g = 9.81; +mu = 0.55; +F = 150; +m = 25; + +% Define the function to be solved for. Is the angle specified in radians or degrees? +%angle=@(-asin(g/(sqrt(mu^2)))+atand(mu)+((4*n2+1)*pi)/2 +func=@(angle) mu*m*g/(cosd(angle)+mu*sind(angle))-F; + +% THINK, what makes range sense for angle? +lower_bound=-90; +upper_bound=90; + +% Plot your function. Does plotting give you an idea about where the root is? +hold on +fplot(func) +fplot(0) + +% Finaly solve for the root using the bisection script given to you, which can be called as: +%[root, fx, ea, iter] = bisect(func, lower_bound, upper_bound); + +angle = bisect(func,lower_bound, upper_bound) + +% These variables will be checked for your final answer: +%angle = % the angle in degree that solves this problem +fx = func(angle) % the function value at your solved angle +ea = ea % the bisection error estimate +iter = iter % the number of bisection iterations
\ No newline at end of file diff --git a/Scripts/hw2/hw2.m b/Scripts/hw2/hw2.m new file mode 100644 index 0000000..6736240 --- /dev/null +++ b/Scripts/hw2/hw2.m @@ -0,0 +1,60 @@ +% Homework:2 +% Author: Christian Kolset +clear + +%% Part 1 + +% Function parameters +q0 = 10; +R = 60; +L = 9; +C = 0.00005; + +% Use linspace to create an array of 100 points between 0 and 0.8 +t = linspace(0,0.8); + +% Calculate the values of q +q = q0.*2.718.^(-R.*t/(2*L)).*cos(sqrt((1/(L*C))-(R/(2*L))^2).*t); + +% Plot q vs t +hold on +subplot(2,1,1) +plot (t,q,'b--*') +title('Capacity vs Time Graph') +xlabel('Time') +ylabel('Charge') +%legend('Charge','Time') +hold off + +% Make the capacitor 10x bigger +q2 = q0.*2.718.^(-R.*t/(2.*L)).*cos(sqrt((1/(L.*10.*C))-(R/(2.*L))^2).*t); + +% Plot q2 vs t +hold on +subplot(2,1,2) +plot(t,q2,'rs:') +title('10x Capacity vs Time Graph') +xlabel('Time') +ylabel('Charge') +%legend('Charge','Time') +hold off + +%% Part 2 +%{ +% Given experimental data +t_exp = 10:10:60; +c_exp = [3.4 2.6 1.6 1.3 1.0 0.5]; + +% Expected function +t_func = 0:0.5:70; +c_func = 4.84*2.718.^(-0.034*t_func); + +% Plot +hold on +plot(t_exp,c_exp,'rd:') +plot(t_func,c_func,'g--') +xlabel('Time [minutes]') +ylabel('Concentration [ppm]') +legend('Experimental','Predicted') +hold off +%}
\ No newline at end of file diff --git a/Scripts/hw3/hw3.m b/Scripts/hw3/hw3.m new file mode 100644 index 0000000..7caf82c --- /dev/null +++ b/Scripts/hw3/hw3.m @@ -0,0 +1,37 @@ +% HW3 +% Author: Christian Kolset +% Date: 20.5.21 + +%% HW3 +% Specify the variables needed to solve this problem (ie. height of each section, diameter, radiaus, ...) +% It is alwasy easier to work with variables (diameter_cyl = 25) than to use numbers everywhere, since a +% diameter indicates something specific but the number 25 could mean anything +diameter_Bot = 25; +r_bot = diameter_Bot/2; +h_cone = h-19; + +% Specify the height of the water +h = 20 +% You can comment / uncomment lines below for testing. This will overwrite the previous line for h = 20. +% For submission, make sure all of the following lines are commented out and h = 20! (OR IT IS MARKED AS WRONG) +%h = 5 +%h = 19 +%h = 47 +%h = -1 + +% Now compute the volume. Using conditional statments you will want to first check the height makes sense, +% and then solve the volume depending on what portion of the tank has been filled. +% Make sure that your volume is stored in the variable v! (OR IT WILL BE MARKED AS WRONG) +% You may find it more convenient to move v around in you code, it is only given here to indicate what variable to use. +%v = + +v_cyl = @(h,r_bot) pi*r_bot^2*h; +v_truncone = @(h_cone,r_bot) 1/3*pi*h_cone*(r_bot^2+r_bot*((h_cone+16.625)/1.33)+((h_cone+16.625)/1.33)^2); + +if h > 19 + v = v_cyl(19,r_bot)+v_truncone(h_cone,r_bot); +else + v = v_cyl(h,r_bot); +end + +fprintf(1,'Volume of Tank: %6.2f\n',v)
\ No newline at end of file diff --git a/Test.txt b/Test.txt deleted file mode 100644 index a9bfdc2..0000000 --- a/Test.txt +++ /dev/null @@ -1 +0,0 @@ -Test text for first commit. |